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ABSTRACT 


The  respome  of  a  hollow  circular  cylindrical  shell  of  arbitrary  thickness,  in 
either  an  elastic  or  a  viscoelastic  medium,  to  transient  dilatatinnal  *ud  shear 
waves  (and  their  superposition)  is  presented.  The  solution  is  valid  within  the 
scope  of  the  linear  theory  of  elasticity  or  viscoelasticity.  The  technique  for 
obtaining  the  solution  relies  upon  1)  the  construction  of  a  train  of  incident  pulses 
i rom.  steady  state  components,  where  each  pulse  represents  the  time  history  of 
the  transient  stress  in  the  incident  wave,  and  2)  the  existence  of  a  physical 
mechanism  that,  between  pulaea,  restores  the  disturbed  particles  of  the  cylinder 
and  the  surrounding  medium  to  an  unstrained  state  of  rest. 


The  influence  on  the  cylinder  response  of  the  following  factors  is  discussed: 
liner  thickness,  cylinder-medium  impedance  mismatch,  viscoelasticity  in  the 
medium,  and  incident  wave  form  ( step  pulse,  rectangular,  triangular,  linear 
rise- exponential  decay). 


u 


CONTENTS 

Section 

Page 

I. 

General 

Discussion  of  Program 

I 

A. 

Program  Logic 

1 

B. 

Input  Specification 

6 

V  * 

Flow  Charts 

.  9 

II. 

Discussion  of  Input  Parameters 
and  Examples  of  Preparation  of 

Input  Data  Cards 

15 

A. 

General  Remarks 

15 

B. 

Calculation  of  Auxiliary 
Parameters 

15 

0. 

Examples 

19 

D. 

Special  Considerations  for 
Early  Time  Response 

34 

III. 

FORTRAN  Listing  of  Program 

35 

Distribution  69 


SECTION  I 


QENERAL  DISCUSSION  OF  PROQRAM 


A.  Program  Logic 

The  program  presents  a  rather  formidable  appearance,  principally  be¬ 
cause  of  it?  length,  but  the  central  idea  Is  quite  simple  -  for  each 
of  the  triples  (*,  •  >x  )  specified  by  the  input,  the  terms  of  a 
double  series  are  computed  and  summed.  Actually  there  are  five  dis¬ 
tinct  series  for  each  (p,  t  ,  x) ,  and  since  each  series  is  a  double 
series,  the  individual  terms  can  be  regarded  as  dependent  on  six 
parameters,  viz,  e  f  9  ,  and  x ,  two  summation  indices,  and  a  final 
Index  to  distinguish  among  the  five  sums.  It  is  this  six-way  index¬ 
ing  that  lies  at  the  root  of  the  difficulties  encountered  in  prepar¬ 
ing  the  program. 

Perhaps  the  most  straightforward  approach  would  have  run  along  the 
following  lines:  choose  p,  9  ,  and  x;  then  generate  the  successive 
terms  of  the  sum  for  p,  9$  and  x  and  perform  the  addition.  Of  course, 
five. sums  must  be  provided  for,  but  this  could  be  done  without  undue 
complication.  The  decisive  consideration  against  such  a  procedure  la 
the  unacceptably  high  amount  of  computer  time  that  would  be  required. 
The  calculation  of  an  individual  term  has  much  in  common  with  that  of 
many  other  terms,  and  it  would  be  pointless  to  ignore  the  possibility 
of  exploiting  the  similarity.  Although  we  have  not  attempted  to  form 
a  close  estimate  of  the  time-saving  made  possible  by  the  procedure 
we  adopted,  a  factor  of  100  appears  to  be  conservative.  Denoting  the 
summation  indices  by  p  and  n,  the  over-all  arrangement  is  as  follows: 
p  is  fixed,  then  p ,  and  then  for  each  n  we  form  the  terms  of  the  five 
sums  associated  with  each  f  and  x  and  sum  over  n;  p  then  tako  on  its 
next  value  and  when  p  has  attained  its  final  value,  p  is  incremented. 
If  k  Is  the  indsx  used  to  distinguish  among  the  five  sums,  then  the 
order  of  variation  of  parameters  ic  p ,  p  ,  n, t  ,  x ,  and  k,  with  the 
understanding  that  p  is  the  "most  permanent "  index  and  k  the  most 
rapidly  changing.  This  arrangement  is  not  necessarily  claimed  to  be 
optimal,  but  it  seemed  to  represent  a  feasible  compromise  between 
the  (sometimes  conflicting)  requirements  of  speed  and  convenient  use 
of  storage. 


Our  discussion  of  the  program  is  confined  to  a  few  major  items  in  the 
belief  that  a  detailed  description  would  be  inordinately  difficult  to 
follow.  As  a  supplement,  we  have  included  two  flow  charts  (at  dif¬ 
ferent  levels  of  detail)  to  facilitate  understanding  of  the  flow  of 
control. 

It  is  widely  recognized  that  the  accurate  summation  of  a  convergent 
series  is  a  severe  problem.  Care  must  be  taken  that  enough  terms  are 
added  and  also  that  the  computed  partial  sums  bear  some  resemblance 
to  the  corresponding  exact  partial  sumB.  In  a  strict  sense,  both 
problems  are  Insoluble  in  a  general  setting;  for  a  computer  program 
can  "see"  only  a  finite  number  of  terms,  and  unless  these  termB  are 
in  some  sense  "typical"  one  cannot  be  certain  that  significant  terms 
may  not  have  been  omitted.  Furthermore,  the  computed  sum  of  even 
three  numbers  may  bear  no  relation  to  the  true  sum,  so  that  even  when 
a  particular  partial  sum  is  itself  adequate,  the  computed  partial  sum 
may  have  no  correct  digits.  These  considerations  indicate  a  need  for 
caution,  and  we  have  acted  accordingly. 

Our  program  faced  two  additional  complications  in  that  the  series  in¬ 
volved  are  double  series,  and  in  the  need  for  "simultaneous"  summa¬ 
tion  of  a  number  of  series.  In  describing  the  process,  we  attempt  to 
maintain  a  degree  of  clarity  by  momentarily  neglecting  the  fact  that 
many  series  are  being  treated  at  once. 

There  are  five  Input  parameters  that  govern  the  summation  of  the 
series.  These  are  NACRCY ,  KAPP,  M*TAP,KAPN,  and  M$TAN.  N AC RCY  stip¬ 
ulates  a  negative  power  of  10  that  is  to  be  an  approximate  bound  of 
the  relative  error  (NAORCY  «  3  indicates,  for  example,  that  an 
attempt  is  to  be  made  to  keep  the  relative  error  below  the  level  of 
10~3,  or  roughly,  to  attain  three  correct  significant  digits).  KAPP 
is  the  maximum  value  that  p  (a  summation  index)  will  be  allowed  (if 
satisfactory  convergence  appears  to  take  place  earlier,  p  may  stop 
short  of  KAPP) .  M^TAP  Is  the  number  of  consecutive  "small"  terms 
required  to  halt  the  addition  ("small"  is  defined  later).  For  each 
p,  we  sum  over  the  index  n,  up  to  a  maximum  of  KAPN  terms,  with  the 
provision  that  if  M^TAN  consecutive  "small"  terms  appear  first,  the 
summation  over  n  is  stopped. 

Associated  with  each  double  series  is  an  Integer  that  plays  a  dual 
role.  Por  present  purposes,  designate  the  series  by  D2  and  the  cor¬ 
responding  integer  by  K2.  The  primary  function  of  K2  is  to  preserve 
a  record  of  the  maximum  size  reached  by  any  of  the  terms  or  partial 
sums  of  D2.  Explicitly,  K2  starts  at  zero,  and  during  the  summation 
it  maintains  a  value  of  128 (m  +  128)  where  2™  is  the  smallest  power 
of  2  which  is  larger  than  the  terms  and  partial  sums  of  D2  so  far  en¬ 
countered.  The  number  128,  which  is  Itself  a  power  of  2,  is 
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convenient  for  use  on  a  binary  machine,  and  the  multiplication  of 
m  +  128  by  128  serves  to  zero  out  the  low-order  7  binary  digits  of 
K 2:  these  digit  onnltlnns  er*e  lined  fnr>  snnthsii  purpose.  A  nearly 
generated  term  is  regarded  as  "small"  if  either  of  two  conditions 
obtains : 

1)  the  term  is  smaller  than  T0L  times  the  corresponding 
partial  sum;  here  T0L  =  .1  x  iq-NACRCY 

2)  the  tern  is  smaller  than  10'?  times  the  power  of  2  (m) 
specified  by  K2.  For  "row"  summation  (i.e.,  over  the 
index  n,  rather  than  p),  the  factor  .5  x  10"?  la  used 
In  place  of  10'^. 

Before  elaborating  on  this  definition  of  "small",  we  note  that  the 
secondary  role  of  K2,  for  which  its  low-order  7  binary  digits  are 
used,  is  to  keep  track  of  the  number  of  consecutive  small  terms  en¬ 
countered;  when  this  integer  reaches  M^TAP  (or  tyiTAN  for  n-summa- 
tion) ,  the  addition  and  generation  of  terms  for  the  sum  is  stopped. 

Roughly  speaking.  Criterion  l)  states  that  the  term  has  no  serious 
influence  on  the  corresponding  partial  sura.  When  this  Is  true  of 
several  consecutive  terms,  it  furnishes  evidence  that  the  remainder 
is  small,  or  that  the  partial  sum  is  adequate;  an  extra  factor  of  .1 
is  Included  so  as  to  be  conservative.  Criterion  2)  is  an  assertion 
that  the  term  is  negligible  compared  to  some  of  the  numbers  to  which 
it  is  added  -  not  necessarily  the  most  recent  partial  sum,  which  may 
(conceivably)  be  substantially  smaller  than  some  that  precede  It. 

The  Intent  or  this  condition  is  to  avoid  the  inclusion  of  many  ad¬ 
ditional  terms  in  a  futile  pursuit  of  better  accuracy;  for  relatively 
low  values  of  NACRCY,  Criterion  l)  will  ordinarily  come  into  play  be¬ 
fore  2),  but  if  NACRCY  is  6  or  7,  the  second  criterion  la  more  real¬ 
istic  and  should  be  dominant. 

When  all  the  summation  is  complete,  these  "descriptive  integers", 
i.e.,  the  K2’s,  are  redefined  so  as  to  estimate  the  relative  accuracy 
obtained.  The  sign  la  positive  if  the  p-serles  convergence  appeared 
satisfactory,  i.e..  If  I^TAP  consecutive  small  terms  were  found; 
otherwise  the  sign  is  negative.  The  magnitude  of  the  new  K2  cannot 
be  greater  than  NACRCY  and  Is  often  less;  a  value  of  m  Indicates  that 
the  relative  error  Is  likely  to  be  less  than  10'ra. 

If  the  final  value  of  a  sum  la  considerably  smaller  than  some  of  the 
intermediate  values,  the  associated  K2  will  be  somewhat  smaller  than 
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would  otherwise  be  the  case.  Cancellation  will  have  taken  place, 
and  this  will  usually  mean  reduced  accuracy,  but  by  "rer.embering"  the 
larger  value,  it  is  possible  to  assess  how  damaging  the  cancellation 
has  been. 

There  are  three  short  subprograms  used  to  monitor  the  convergence 
oroeedure.  IBI5  lias  two  arguments,  A  and  N,  and  its  purpose  is  a 
magnitude  comparison.  N  is  thought  of  as  a  "descriptive  inr-eger"  so 
that  its  high  part  represents  a  power  of  2.  If  this  power  cf  2  is 
as  large  as  the  magnitude  of  A,  then  IJ3IG(A,N)  is  given  the  value  N; 
otherwise  M  is  determined  to  satisfy  2K_2-  -  IAI  <  and  the  high  part 
of  IBIG(A,n)  is  set  to  128(M  +  128)  with  the  lew  part. (7  least,  sig¬ 
nificant  bits)  taken  from  K.  TW^K  has  the  single  argument'  K,  and 
Tty0K(K)  is  the  floating-point .number  2^' if  the  high  part  of  K  is 
128(M  +  128) .  LL0W  aleo  hab  the  single  argument  K,  and  LL^W(K)  is 
Jus1.;'  the  integer  given  by  the  low  >rdn*r  7  bits  of  X.  There  subpro- 
grams  are  given  in  .F^PTRAN,  but  equivalent  unsembL;-  language  programs 
are  easily  written  and  operate  much  more  quickly. 

Aside  from  this  convergence  procedure,  most  of  the  coding  is  fairly 
straightforward  and  not  strikingly  novel.  Perhaps  the  most  unusual 
feature  Is  the  Bessel  function  computation,  and  even  here  the  method 
has  become  almost  standard.  The  program  requires  Jq(z),  J-^(z),. 
for  a  given  (possibly  complex)  argument  z;  there  1b  also. a  need  for 
Yq(z)  , 'YAz) , ,  and  for.  some  scaled  Bessel  functions  that  we  refer 
to  as  and  Y^,  but  each  Y  sequence  is  easily  obtained  from  the  cor¬ 
responding  J  sequence  and  the  a^'s  are  computed  in  virtually  the  same 
manner  as  the  Jn'e,  so  we  discuss  only  the  computation  of  Jq,  Jj 
etc.,  for  a  given  z.  The  basic  tool  used  is  the  recurrence  formula, 

•W*>  -T 

For  a  fixed  z,  this  formula  can,  in  principle,  be  used  to  generate 
Jn+^(z)  for  n>l,  if  Jq(ss)  and  J^z)  are  known.  Put  the  J  sequence 
decreases  very  rapidly  with  increasing  n,  while  the  complementary  Y 
sequence  increases  very  rapidly.  Unless  the  values  taken  for  JQ  and 
are  exact  and  all  subsequent  calculation  is  also  exact,  the  com¬ 
puted  J  sequence  will  contain  an  initially  small,  but  eventually  dom¬ 
inant  component  of  the  Y  sequence.  For  typical  values  of  z  and  good 
single-precision  approximations  of  J Q  and  J^,  we  often  find  a  computed 
J10  t0  be  worthless. 

It  is  possible  however,  to  turn  this  behavior  to  advantage;  if  Jq, 
Jl,....Tn  are  desired,  one  selects  a  suitable  m>n  and  defines  »  0, 
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pm-l  =  *  with  *  small  but  arbitiary.  The  recurrence  formula  is  used 
in  reverse  order  to  generate  the  sequence  Fm_2>  p..  01  po* 

The  F  sequence  then  satisfies  F^  =  bJ^  +  0Y^  for*  ome  a  ard  0  inde¬ 
pendent  of  kj  in  particular,  up  v-aw  Fm  -  0  -rJm  :  07^,  so  that 

JL  =  and  we  can  then  write  a  with  relative  error  satis- 

a  Ym 

fying 


If  m  is  considerably  larger  than  k,  we  will  have|jra  «|  |and  also 
lyk!<<;|  so  that  the  relative  error  is  the  produce  of  two  small 
factors,  and  it  is  reasonable  to  assume  P^  “  at  .least  for  com¬ 
putational  purposes.  In  practice  it  is  sufficient  to  take  m  equal  to 
n+30  to  produce  good  single  precision  accuracy  for  fairly  small  z, 
say  |z|<  50,  To  complete  calculation  of  the  J's,  we  need  only  de¬ 
termine  the  value  of  a  .  This  is  done  by  using  the  relation: 

wp(z)  +  2Jg(z)  +  2J^(z)  +  •<••••  *»  1. 
in  the  form 


P0(s)  +  2Pg(z)  +  2F^z)  +  . -  («)(D- ■«  • 

If  convergence  of  this  series  is  overly  slow  -  and  this  happens  very 
rarely,  if  at  all  -  then  we  substitute  a  direct  calculation  of  Jc{z) 

bo  yield  «  This  algorithm  is  somewhat  prone  to  overflow,  so 

Jo 

provision  is  made  for  rescaling  the  P  sequence  if  necessary.  If  z  is 
complex,  then  *  is  ordinarily  complex  and  caution  must  be  exercized 
ir  dividing  by  «  .  In  the  event  a  («,  a  +  ib)  is.  large,  one  should 
avoid  forming  a2  +  b  ,  which  may  easily  overflow: when  a  and  b  do  not. 
Recurrence  techniques  for  generating  Bessel  functions  have  been  dis¬ 
cussed  by.  Goldstein  and  Thaler.  * 

Another  problem  of  computational  interest  is  the  solution  of  a  set  of 
linear  equations.  The  subroutine  LINjtfX#  has  been  used  to  accomplish 
this  task.  It  employs  a  direct  Gaussian  elimination  scheme,  with  so- 
called  partial  pivoting.  This  means  that  at  the  kth  stage,  the  pivot 


*  Goldstein,  M.  and  Thaler,  R.  M.  "Recurrence  Techniques  for  the 
Calculation  of  Bessel  Functions",  Mathematical  Tables  and  Other 
Aids  to  Computation,  Vol.  13,  65-b8,  1959,  PP-  102-108. 
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is  cho«en  as  the  largest  element  in  the  k^h  column,  except  that 
above-diagonal  elements  are  not  eligible  for  consideration.  An  in¬ 
terchange  ox*  rows  tnen  takes  place,  if  necessary,  to  place  the  pivot 
in  the  diagonal  position  of  column  k.  We  then  subtract  multiples  of 
row  k  from  the  subsequent  rows,  so  as  to  introduce  zeros  into  the 
lower  part  of  column  k.  After  n-1  steps  of  the  algorithm,  the  origi¬ 
nal  system  of  equations  has  been  replaced  by  an  upper  triangular 
system  whose  solution  presents  no  special  difficulties. 

The  roles  of  the  subprograms  IBIQ,  TW.0K,  LI^W,  and  LIK0Y0  have 
already  been  described.  For  convenience,  we  list  the  remaining  sub¬ 
routines  and  briefly  indicate  their  functions: 

BESSEL  calculates  the  sequence  Jo(z),  Ji(z),  «?n(2) 
for  real  or  complex  z. 


TILDE 


£omputes  tjje  scaled  Bessel  functions 
J1(t) . Jn(z);  and  also  Y0(*>»  M.(*) 


J0(*) 


These  functions  are  defined  In  Vol*  I,  Appendix  A, 


equations  (A. 13). 


JTIL0  is  used  by  BESSEL  AND  TILDE  to  supply  a  value  for  J0(z). 

ZRHSP  and  RHSP  generate  the  right  sides,  M„  of  the  systems  of  ■ 
equations  for  n  =  0  and  n  >0,  respectively,  if 
the  wave  is  dilatational. 

ZRHSM  AND  RHSM  generate  these  right  sides  for  the  shear  wave. 

Aside  from  these  subroutines,  the  program  employs  only  the  standard 
library  tape  subprograms. 

B.  Input  Specification 

The  first  card  has  the  format  (A6,  1713).  It  should  contain  a  six- 
character  code  to  identify  the  run,  followed  by'  13  integers.  The 
code  may  include  any  combination  of  the  10  decimal  digits,  the  2 6 
letters,  and  the  various  special  characters,  including  blanks.  It  is 
reproduced  on  the  output  but  has  no  other  effect ,  The  13  integers 
are,  in  order,  K0DEN,  KTYPE,  K0D0P,  KC0,  MAX0UT,  KAPP,  MgTAP,  KAPN, 
M#TAN,  NACRCY,  NR0  ,  NTH,  and  NMAJ;  they  should  appear  in  columns 
7-9,  10-12,  13-15,  etc.  The  significance  of  KAPP,  MfifTAP,  KAPN,  M#TAN 
and  NACRCY  has  already  been  indicated.  K0UEN  is  taken  as  0,  1,  -1, 


accordingly,  as  the  run  is  to  describe  a  Bimple  dilatational  wave,  a 
simple  shear  wave,  or  a  combined  wave,  respectively.  _  The  value 
KTYPE  =  1  is  used  for  a  rectangular  wave,  KTYPE  =  2  for  a  triangular 
wave,  and  ICTYTE  -  5  an  exponential  aecay  wave.  For  a  combined 
wave  (K0DEN  =  -l),  this  value  of  KTYPE  applies  only  to  the  dilata¬ 
tional  part  and  «  new  value  is  to  be  furnished  for  the  shear  wave; 
we  shall  return  to  this  consideration  later.  K0D0?  should  be  given 
as  1  if  a  Lanczos  factor  (/tee  Vol.  I,  Section  VT.  2)  is  to  be  in¬ 
cluded,  and  as  0  otherwise.  With  KO0>O,  there  is  a  large  volume  of 
intermediate  output  (see  the  sample  output).  If  KC0  <  0,  the  inter¬ 
mediate  output  is  suppressed.  Similarly,  if  MAX^UT  /  0,  the  maxima 
(of  the  stress  quantities)  with  respect  to  P,  t,  and  X  and  the  over¬ 
all  maxima  are  computed  and  printed.  If  MAX0UT  *  0,  this  output  is 
not  given.  Ni^0  is  the  number  of  ?  values  to  be  used.  NTH  i‘;  the 
number  of  G  values  explicitly  given.  NMAJ  is  the  number  of  'major" 
x  values  given,  but  it  is  possible  to  insert  equally  spaced  sub¬ 
division  points  between  adjacent  X  values,  as  discussed  below,  in  con¬ 
nection  with  the  array  INTERX . 

It  should  be  noted  that  all  the  floating-point  input  is  to  be  given 
in  (E14.7)  format.  This  applies  in  particular  to  the  paramete/s  E., 
Eg,  ...E  ,  and  to  the  values  of  *  ,  the  NTH  values  of  G,  an!  trie 

NMAJ  values  of  x,  which  follow  the  fir*t  card  in  the  order  given. 

If  NMAJ  ■=  1,  the  next  card  is  omitted,  but  otherwise  there  Bhould 
then  ba  a  card  containing  NMAJ-1  integers  in  (2013)  format.  These 
integers,  called  INTERX,  indicate  the  number  of  intermediate  X  pointB 
to  be  inserted  between  the  given,  so-called  "major"  X  values.  Speci¬ 
fically,  INTERX (i)  stipulates  the  number  of  new  points  between  XMAJ(I) 
and  XMAJ(I+1).  ,  v 

If  KTYPE  -  1,  this  completes  the  input  for  a  simple  dilatational  or  ' 
shear  wave,  or  for  the  dilatational  part  of  a  combined  wav  >.  If 
KTYPE  ■  2  or  3»  two  more  (E14.7)  cards  are  required,  wi-h  values  of 
CTBIN  and  CAYIN  or  EMIN,  respectively.  No  additional  input  is  needed 
for  either  of  the  simple  waves,  but  If  a  combined  wave  is  to  be  de¬ 
scribed,  several  cards  follow  to  establish  the  desired  relation  be¬ 
tween  the  two  parts.  First  is  an  Integer  card  containing  LEADX, 

KEYX,  and  NTH  values  KEYTH.  LEADX  specifies  the  XMAJ  value  that  is 
to  become  the  first  it  for  which  the  shear  wave  calculation  is  to  be 
performed,  i . e. ,  XMAJ( LEADX)  serves  as  this  x.  Similarly,  XMAJ (KEYX) 
designates  X*,  the  nondimens ional  time  delay  between  pulses.  Next  Is 
an  integer  card  containing  the  values  of  KTYPE  and  KCJ0  that  are  to 
apply  to  the  shear  wave  calculation,  and  then  three  cards  containing 
PHI,  AMPFAC,  and  a  new  Ell  value.  If  KTYPE  -  2  or  3,  two  additional 
cards  are  needed,  with  the  new  valves  of  CTBIN  and  CAYIN  or  EMIN, 
respectively.  * 


■P'; 


<^:,r  • 
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The  correspondence  between  FORTRAN  input  symbols  and  the  actual  para* 

Uocu  xil  1 1C  aiiaiy  SXD  1  U1  ,LOWB  I 


E(l)  -  h/R 
E(2)  -  -rA 


— V 

w 


E  /E 
m 


E(6)  -  ^ 

E(7)  “  b02 
cdm 

E(8)  -  0X 

E(9)  -  fyn 
E(l°)  -  r  S 


E(ll)  ■  cdmto 
b 

PHI  -  $ 


AMPHAC  - 


T.iner  thickness  tc  mean  radius  ratio 

Inverse  of  non-dimensional  1/2 
period  of  the  wave  form 

Poisson's  ratio  for  the  liner 

Poseion's  ratio  for  the  medium 

Medium  to  liner  ratio  of  Young’s 
Modulii 

r 

Liner  to  medium  ratio  of  mass  density 


Inverse  of  non-dimensional  relaxation 
time  for  shear  stress 

Stress  relaxation  time  versus  uniaxial 
strain  relaxation  time 


Stress  relaxation  time  versus  strain  re¬ 
covery  time  for  uniaxial  strain 

l-  •  '  ”  - .  1  V. 

Shear  stress  relaxation  time  versus  shear 
strain  relaxation  time  ......  .. 


Nondim^nBional  rest. time 


Angle  between  incoming  dilatational  and 
shear  wave 

Amplitude  ratio,  shear  wave  to 
dilatational  wave 


':  *  ■  ;  ,  ■ »  ■’  *  ( '  -  £  ‘  '  ^  ■ 

1 

k 

•  ■'  v  -  '<  ■ .  .  . 

“  f:  \T-‘vri  .—r-—  -  ■  -- 


CAYIN  =  k 


CTBIN  - 

b 


EMIN  -  k 


Integer,  describing  the  decay  time  of 
triangular  wave  form;  decay  time  = 

k(s(yp) 

Non-dimensional  rise  time  for 
triangular  or  linear  rise-exponential 
decay  wave  form 

1cdm  dilatational  wave 
°tm  shear  wave 

Parameter  associated  with  linear  rise  - 
exponential  decay  waveform  and  whioh 
defines  the  time  at  which  the  stress  in 
the  exponential  decay  portion  of  the 
pressure-time  history  is  negligible, 

e'k  <<  1. 


C.  Flow  Charts 

This  subsection  contains  two  flow  charts,  intended  to  assist  in 
the  understanding  of  the  program's  flow  of  control.  Chart  I  is 
relatively  brief;  It  attempts  to  convey  a  sense  of  the  over-all 
structure  of  the  procedure  and  keeps  the  use  of  special  symbols 
to  a  minimum.  More  details  are  presented  in  Chart  II. 
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Compute  Print  Maxima 


11 


Adju»t  Parameter*  and  Set 
KEY  *  -I,  JUMP  -  0 


CHAU  II  (C-  TIHUtt) 
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SECTION  II 


DISCUSSION  OF  INPUT  PARAMETERS  AND 
EXAMPLES  OF  PREPARATION  OF  INPUT  DATACARDS 


A.  General  Retnarka 

The  efficient  use  of  the  computer  program  depends  in  part  on  the  proper 

choice  of  some  key  parameters.  These  are  KAPP,  KAPN,  MJITAP,  MJPfAN, 
NACRCY,  Cfflto/b,  and  b/Cfap.  In  all  the  data  presented  in  Vol  I, 

MJiTAP,  M^rAN,  NACRCY  were  taken  as  4,  6,  4,  respectively.  These  values 
were  considered  to  he  adequate  for  obtaining  the  plotting  accuracy  de¬ 
sired.  Although  increasing  the  value  of  these  parameters  should  in¬ 
crease  accuracy  (and  also  increase  the  computer  time),  the  exact  lim¬ 
itations  and  effect  of  such  an  Increase  has  not  been  studied. 

The  accuracy  of  any  computation  is  also  dependent  on  the  choice  of 
KAPP  since  KAPP  controls  the  Fourier  representation  of  the  wave  form. 

In  general,  KAPP  wab*.  chosen  qa  101  for  the  rectangular  incident  wave 
form,  200  for  the  triangular  wave  form  and  greater  than  300  for  the 
linear  rise  exponential  decay  wave  form.  These  values  were  adequate, 
in  moqt  cases,  for  obtaining  approximately  four  figure  accuracy  for 
all  the  quantities  of  interest.  KAPN  was  chosen  as  50  for  all  the  re¬ 
sults  presented  In  Vol  I.  However,  this  value  is  conservative  and  in 
most  cases  it  can  be  decreased  considerably.  The  determination  of 
b/cdoT  and  c^to/b  is  discussed  in  the  next  subsection. 

B.  Calculation  of  Auxiliary  Parameters 

The  significance  of  the  parameters  controlling  "rest  time"  and  the 
wave  form  were  discussed  in  Sec.  V.l  and  Appendix  B  of  Vol.  I.  In  the 
present  subsection,  the  procedure  for  determining  the  value  of  these 
parameters  is  discussed  for  each  wave  form.  It  is  assumed  that  mater¬ 
ial  and  geometric  properties  of  the  liner  and  medium  are  given  (l.e,, 
cdm*  cta>  ?>)  and  that  only  non-dimensional  parameters  are  discussed. 

To  facilitate  the  discussion,  components  of  Fig. '2b,  Vol  I,  showing 
the  various  wave  forms  and  associated  parameters ,  are  reproduced  im¬ 
mediately  prior  td  the  discussion  of  the  wave  form  being  considered. 
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mctanoulaa  wive  form 

Til 


cdm  *er  incoming  dllatorioncl  nova 
elni  I  of  ‘neemino  tutor  watt. 


As  mentioned  in  Appendix  B,  Vol  I,  any  two  of  the  parameters  b/crimT„: 
Cmtc/b,  c_  A/b  may  be  chosen  arbitrarily.  In  preparing  input  data  for 
the  reported  problems  it  was  most  convenient  to  choose  cmt0/b  and 
cm  A/b,  the  dimensionless  "rest  time"  and  pulse  length,  respectively, 
and  then  solve  for  b/cqmT  from 


b 

cdm^ 


cnv/'Qdm 


2( 


<^0 


>  +  ( 


(1) 


which  is  obtained  from  Eq.  (5.s)  Appendix  3,  Vol  I. 

To  demonstrate  this  computation  for  a  dilatatlonal  Incident  wave,  we 
choose 


cdmfco 

b 


20 


■  .  20, 

b 

and  find  from  Eq.  1,  that  b/c^T  0.01666667. 


For  most  of  the  problems  that  were  studied,  2cmtq/b  *40  (with  20  £ 
CdmTA>  £  200)  was  sufficient  to  ensure  rest  conditions  before  the 
pulse  strikes  the  cylinder. 


-rr^m 


TFiANGOLAR  WAV'  FORM 
Fm  Ut 

^.1 _ 

4 


0  I - 1  ■  I,  _ i_| 

k—  *o  *i  *  .-  •»  >U|  •  i0h. ..— u.  ij  J 


e  .  iC,)m  (O'  Incoming  dllotitlonol  «o«r 
41  letm  (or  incoming  thtor  wovt 


The  parameters  that  define  the  wave  form  and  rest  time  are 
k>  b/cdmT-'  and  cmtoA  Al  discussed  in  Appendix  B,  Vol.  I,  only  three 
of  these  parameters  can  be  chosen  arbitrarily.  Since  it  Is  natural 
that  cm(ti-t0)/b  and  k  are  prescribed,  we  are  free  to  choose  either 

a)  b/cdmT'  and  calculate  c^o/b,  or 

b.)  cmtoA»  and  calculate  b/c^T. 


Again,  for  the  triangular  wa\ i  form  it  was  most  convenient  to  choose 
cratc/b  along  with  k  and  cm(t1~tc)A»  and  find  b/c^T  from 

b  cm/cdm 

cdmT  "  2cmfc0  +  (k  +  1)  cm(t:rtc)  ' 
b 

where  Eq.  2  la  obtained  from  Eq,  (B.6)  Appendix  B,  Vol  I. 

Thus,  for  a  dllatational  incident  wave,  given 


cm(tl~to)/b  ’  1 


we  choose 


rise  time 

and  decay  time  of  pulse, 
2cmto  .  40, 


and  i-nd  from  Eq.  2  b 

cdm^ 


0.021739. 
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LINEAR  RISE-EXPONENTIAL  DECAY  WAVE  FORM 


tim  for  incoming  dilototionol  wovo 
e  in  lot  incoming  th»or  wove 


Note  *  it  toltcttd  90  that  the  Incident  droit  I*  ntgllgoably  (moll  ft  t*T 


The  parameters  defining  the  wave  form  and  rest  time  are  b/c^mT, 
cmt0/b,  cm(k/«)/b,  cm(ti-t0)/b.  Again,  we  are  free  to  select  only 
three  of  these  parameters  arbitrarily.  Since  cm(k/a)/b  and  cm(tx**t0)/b 
are  always  given,  we  can  select  either  b/cfaT  or  cmt0/b  arbitrarily. 


As  for  the  other  two  wave  forms,  it  is  again  most  convenient  to  choose 
cmt0/b  and  calculate  b/efafl.  Thus,  from  this  wave  form, 


cdmT 


tn 


'dm 


Vo 


vw 


Hence,  for  an  incident  dilatational  wave,  given 


cdo(l/a> 


and 


c  .  (t, -t  ) 
dmv  1  o 


c  .  t 

we  take  k  •  12  and  — r—  ■  30 


and  calculate  b/c^T  ■  0.010989. 


A 


f»' 

% 


.  ir  :  *. 
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C.  Example* 

Example  1  -  Preparation  of  Input  for  a  dilatational  rectangular  wave 
followed  by  a  rectangular  shear  wave. 

Given  geometric  and  material  properties  are: 

Slow  Granite  Medium  (Elastic)  Concrete  Liner 


Em  .  1.0  x  10°  pel 
m 


m 


.25 

y  ■  5.2  slug/ft3 


m 

*»1 


1.0 


•  0.0 

cdm 


Calculated  quantities  are: 

h 


E  -2.5  x  10°  psi 

*  -  .2 

y  -4.5  slug/ft3 
Thickness  h  •  2.391  ft.. 

Mean  Radius  R* 23.81  ft, 

v  r 

or  ■ 

Outer  Radius  b  •  25  ft. 
Inner  Radius  a  •  22.62  ft. 


•  .1 


'£bl.. 4 
.  E 


■=—•  .86538462 


dm 


°tm 


- 


d  ■  vm> 


^  d  - 


Vij)  d  +  vra) 


2d  +  »„)  I'm 


•  5765  ft/sec . 


•  3328  *'V  3  '  . 
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■:V  r_\ 


*VV 


and  we  calculate 

.  -  1  ’  '  1 

.  Jl/  °tm\  /  1  -  °dro  Ag\.  it  lias. 

*  a\W  \b/v  ~  1 

For  an  air  Induced  loading,  traveling  along  the  ground  at  the  super- 
seismic  velocity, U  -8500  ft/sec,  the  angle  between  the  dilatational 
and  shear  wive  is  found  to  be  e  •  19.7°.  If  it  is  aLjumed  that  each 
wave  travels  Approximately  150  ft  before  it  strikes  the  cylinder, 
then  the  shear  wave  strikes  the  cylinder  19  mil  sec  after  the  dilata¬ 
tional  wave.  The  non-dimensional  time  delay  is  therefore  found  to  be 

i:  1 

1 

°dra  4t  4.381. 

— ...  m  m  . 


•Q 

On  empirioal  grounds,  the  amplitude  ratio;  *0  ,  of  the  two  incident 
waves  is  te’ren  to  be  1/3.  To  complete  the  input  information,  p  ,  $  , 
x  points  are  chosen. 

For  the  demonstration  problem,  five  ©  points  were  chosen  for  the 
dilatational  part.  Two  extra  ©  points  are  selected  for  the  shear  and 
combination  part  [specified  by  KEYTH(I)J.  Five  major  stations  for  x 
were  chosen;  one  point  (-1.0)  before  the  wave  strikes  the  cylinder,  to 
check  that  all  stress  and  displacement  quantities  are  approximately 
zero,  and  a  point,  X  •  4.381,  corresponding  to  the  time  the  shear  wave 
strikes  the  cylinder.  It  should  be  mentioned  that  this  latter  point 
is  required  information  for  all  combination  runs.  An  additional  24 
points  are  also  Included  f controlled  by  INTERX(l)].  Only  one  a  point 
ia  used  in  this  case,  p  ■  jo. 9047619,  the  inner  surraoe. 

A  chart  with  all  input  data  for  this  example,  la  shown  on  the  next  two 
pages,  .  i 
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RUN  NO.  I  nOO£N  1  KTTPE  I  KOOSP  I  KCO  JS»*  OUfJ  K  APP |  M  ST**  [  K  APN  [  MS  TAN  INACRCy!  NRO  j  NTH 


1 


'  i 


i 


!  • 
■  ’.'WV, 


i 

1 


* 


Example  2  -  Preparation  of  input  for  an  incident  rectangular  shear 
wave. 

The  preparation  of  input  for  this  case  follows  from  the  discussion  of 
a  dilatational  rectangular  wave  (dilatational  case.  Example  l)  with 
the  following  changes. 


:t-V* 


1)  KODEN  -1 

2)  cdmto 


3) 


ctmto 


ctn/cdm 


cdmT 


(ctm^o  \  +  /  ctmdk\ 


Accordingly,  choosing 

% 

*ti 

T 


♦jit.  <•*  »- 

*  Vi . 


-'St  ‘ 


ctmto 


.  ».  "lth  iad. .  io. 


*f  vs 


■  4V .1  «-■  ■«_  •  -  rrrr  - 


we  find 


n  'V 


cdmT 


•  .01155 


A  chart,  with  all  input  data  for  this  problem  is  on  the  next  page,. 


■  jjj 
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I 


,v  -f  ■ 

.  V.  j 


!  '  l:  '  ‘P 


-  Preparation  of  lnp'it  for  an  incident  triangular  shear 

wave  iMiro  -  2) . 


Medium:  Porous  Sandstone  (Elastic) 


E  -  3.0  x  lO  psl 

m 


p  ' m  0.2 

HI 

7ro  •  t.Wi  slug/ft3 


m 


1.2 


4-  m  i  ,0135-35 

y  ■■  v  ;  - 

n 

Parameters  describing  tne  wave  form: 

'  ‘  -  : 

- Mat*, 


Liner: Concrete 
E  *2.5  x  10gpsl 

*  *0.2 

y  -4.5  slu^/ft3 

I*0-1 

t*  25  ft. 


Rise  time 
Decay  time 


10  mil  sec 
50  Sill  seo 


form  the  parameters  describing  the  triangular  wave 

‘form.  aye. 


.-as  v-  a 


•iV  *'V- 


•  ‘tm'V'o1  A  *  £-5'*68 

k  •  5  * 


Choosing  the  peat  time  as 

t  ’  r.  »•  ni»  u*  ^  n*.  V  . 

'  °t»V'b  *  *> 


««■  . 
v.'t'fCT! 

, V,| '  M" 

Wt.-'ii-'"]  'v  ■  s'  1  ■■ 

■  - ' l>'  ;  •  we  calculate  -"b/c^T  from 


A  l. 

■--v.  *  •* 


>--:1  K . 

. 

ray  | 

fc-f 
il'V' 

te? 

r  %ii  ■ . 
1:^1  - 


"dm 


t 

5S-£ 


itw/°dm 


^  +  ^jotm(tl~^o) 
b 


.0081389. 


,VTo  complete  the  information*  we  choose  the  #  *  #  *  x  points  at  which 
calculations  are  desired.  In  this  example*  three  ©  points*  one  * 
point  and  two  x  major  points  are  chosen. 

■  •  ."  ,o,  *  ■*  , 

The  Input  for  a  dilatatlonal  incident  triangular  wave  form  follows 
from  tjilp  dlscxiaalon  with  the  Indicated  changes: 
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'.«•<  nr*  —  *'•»«  TV  ’I'l'***!?’*'’  1  ir-  - 


'??TT2'&» 


V*  ™ 


1 )  KpDEN  -  0 


°tmVD 


v</> 


II 


3)  VctoT  - 


-dB>.°-  +  (k  + 1)  °*>^.rt°) 


The  following  chart  contains  the  input  data  for  thia  problem: 


sr 


'iSjaKT-- 


RUN  NO  I  KOOEN  KTTP6  KOOSP  KCO  MAX  C'JTl  K  *>>P  [mSTAP|  K APN  I  M STAN InACRCvI  NRO  [  NTH 


Example  4  -  Preparation  of  input  for  an  incident  dilatational  linear- 
ri s e- exponen t ia  1  decay  wave  form  (KTYPE-3). 

The  geometric  and  material  properties  are  the  same  as  Example  3. 
Pai'aiiicl/cro  uceuribing  the  wave  form  are: 

Rise  time  ■■  10  mil  sec 


Therefore, 


and  choosing 


decay  time  1  cr.  .  mAn 

constant  “  a  50  mil  sec 


WvV/*  “  4*1554» 

cdm^1/ *  >  /  b  *  20-78, 


Wo  /  b  '  30' 


cdmT  cdmfco  +  °dm^tl~t0^  °dm*k^tt) 


.0041326. 


It  should  be  noted  that  since  this  wave  fox*m  is  the  most  difficult 
to  represent,  KAPP  must  be  taken  greater,  than  301  in  order  to  ensure 
proper  convergence. 

The  incident  shear  wave  case  is  handled  in  an  analogous  manner,  with 
the  following  exceptions: 


1)  K0DEN  •  1 

2>  V»A  — 


0tmtc/b 

°tm(Vto)/b 


3)  *>/<>*,*  is  calculated  from 


8/  CjjjjT 


ctn/cdm 

•tabO+-et«(Vto>  4  0t»(k/«  1 


The  input  data  for  this  problem  are  given  in  the  chart  on  the  follow¬ 
ing  page. 
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RUN  NO.  |  KOOCN  I  KTYPE  I  KOOSP  I  KCO  TldAX  OUTTkAI**  I*  STAR  I  K  APN  ]  MSTAN  INACNCY]  MftO  I  NTH 


FOR  ATYPt  •  *,S 


* 


SAamuic  5  -  Preparation  or  input  ror  a  rectangular  wave  when  the 
medium  Is  viscoelastic. 


The  preparation  of  Input  for  a  viscoelastic 
from  the  elastic  cases  previously  discussed. 


“i 

°2 


’1 

'TTT' 


medium  follows  formally 
That  is,  while  °°2» 

cdm 


must  be  chosen  appropriately,  we  choose  the  other  input  as previously 
discussed  (using  the  same  values  for  the  other  parameters,  the  pre¬ 
vious  examples  in  which 


cdm  ^2 
b 


1.0 


describe  the  elastic  unrelaxed  case). 

Partly  because  viscoelasticity  slows  down  the  process  of  strain  re¬ 
covery  in  the  particular  medium  chosen,  it  is  much  more  difficult  to 
pick  rest  times  that  are  adequate  to  ensure  that  the  initial  condl-: 
tione  are  satisfied  at  the  time  of  arrival  of  the  pulse.  However, 
this  difficulty  can  be  surmounted  in  most  cases  by  conatructing  each 
viscoelastic  problem  from  two  or  more  individual,  oomputer  runs.  For 
the  first  run  the  duration  of  the  pulse  is  made  small,  thus  increas¬ 
ing  the  allowable  rest  time  and  also  decreasing  the  amount  of  creep 
deformation  in  the  medium.  A  subsequent  run  can  be  one  in  which  the 
desired  pulse  length  or  any  intermediate  length  is  used.  For  each 
subsequent  run  the  initial  conditions  may  not  be  satisfied,  but  at 
later  times  the  results  should  be  acourate  enough  to  match  those  of 
the  previous  run.  Ultimately,  the  response  over  the  desired  pulse 
duration  can  be  constructed  in  this  way. 


An  example  set  of  data  for  a  case  for  which  this  procedure  was  used 
1b  provided  at  the  end  of  thiB  discussion.  It  should  be  noted  that 
the  effects  of  viscoelasticity  on  the  satisfaction  of  initial  condit¬ 
ions  as  discussed  above  were  not  severe  for  the  triangular  and  linear 
rise  exponential  decay  wave  forms  that  were  studied. 

The  procedure  outlined  above  does  not  apply  only  to  viscoelastic 
cases.  It  should  also  be  used  for  an  elastic  case  whenever  difficulty 
is  encountered  in  attaining  sufficient  rest  time  and  when  more  refined 
results  are  desired  at  early  transit  times. 

Example  input  charts  for  a  step  pulse,  in  a  viscoelastic  medium  in¬ 
volving  two  computer  runs,  are  given  on  tho  foilwing  pages: 
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examu  5  nror  data  (cost) 

nun  no.  I  koocm  1  urvPC  I  koosp  I  kco  (maxout{k«pp  1ms  tap  I  k  apn  [  mstan  InacrcyI  mmo 


ADDITIONAL  INPUT 
FON  KTVPC  •  *.5 


I] 


.  Special  Considerations  for  Early  Time  Response  (for  Step  or 
RAnttflMtMil  ftr*  Pul  non  1 

i  .M.  - -  ...  ■■  t 

The  accuracy  of  the  early  time  responses  depends  on  the  choice  of 
KAFF .  Omto/b*  and  b/c<jmT.  The  parameter  cmt0/D  controls  the  "rest 
time"  wnlle  KAPP  affects  the  Fourier  representation  of  the  incident 
wave.  The  parameter  b/c<jmT  affects  the  early  time  response  due  to  a 
"scaling  effect"  that  can  be  described  as  follows:  For  a  finite  num¬ 
ber  of  terms  the  Fourier  expansion  approximates  the  step  discontinuity 
in  the  incident  pulse  by  a  steep  ramp.  The  parameter  hfcdtpT  governs 
the  slope  of  this  ramp  relative  to  the  size  of  the  cylinder.  That  is, 
in  the  scale  of  the  cylinder,  a  decrease  in  the  value  of  this  param¬ 
eter  corresponds  to  a  stretching  of  the  ramp  region  in  the  pulse.  Ac¬ 
cordingly  the  inaccuracies  in  the  cylinder  response  due  to  the  passage 
of  this  ramp  over  the  cylinder  will  occur  over  a  larger  interval  of 
time.  By  a  reverse  process,  the  accuracy  of  the  early  time  response 
can  be  improved,  provided  that  the  increase  in  b/c<amT  does  not  de¬ 
crease  the  value  of  Cnjto/b  and,  hence,  the  rest  time,  below  an  accept¬ 
able  limit.  To  keep  the  rest  time  within  this  limit  may  then  require 
getting  a  refined  early  time  response  with  a  pulse  of  short  duration 
and  then  proceeding  as  outlined  previously,  in  Example  5  to  get  the 
complete  response.  In  most  practical  applications  where  the  maximum 
response  is  of  greatest  Importance,  the  slight  inaccuracies  introduced 
at  early  time  are  acceptable,  and, as  in  the  example  problems  presented, 
the  above  procedure  need  not  be  applied. 
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FORTRAN  LISTING  OF  PROGRAM 


C  MIN  RR09R4N 

C  DYNAMIC  STRESS'S  O'  THICK  CVLINOtRS 

DIMENSION  KCON1I 10, 1901, 01(10,1901, KCON2(9, 10. 1901,02(9*10, 1901 
DIMENSION  LCONII  10i90l  ,  MCONI  <  10,301 ,  NCONII  10,901.  ICCNK  10, 301 , 
1JCON11 10, SOI  •  LC0N213, 10,301 ,  MCONI  (3.10,30 1,  MCONI  1 3, 

3 1C  ONI  I  Si  10.301 ,  AK10.30I,  01(10,30),  CI410.I0I.  91110,901. 

SHI ( 10,301 .  A2I3 ,  10.301 ,  9213, 10,301,  CI(3,10,30I,  02(3*10,301, 
4M2I3,  10,301,  JCON213, 10,301 
DIMENSION  f  (3.101, J'(3,tOI,  COSRTI 30 1  ,Sf  NRT(  301 

DIMENSION  INAJI10I  •  KEY  THI 10 1 ,  INTIKI  19 1 
DIMENSION  TADI30I , TSII30I ,TCC 1 10 1 ,TCC( 301 , THTI 30 1 - 
OIM'NSION  RH0I3I .TMETAI10I ,XBAR(30I ,T I TLE t *1 I , I ( 111 
DIMENSION  REJI100I,  ENJI100I,  MASTEIIOOI,  EiOlRIlOOl,  BJQIB(IOO), 
1BJ01I100I,  B  J024  1001 ,  9RJ3I100I  •  9IJ3I100I.  BMJ4IIOOI,  0U4I100I, 
I0JQ1RI 1001 ,  BJB2RI 1001  •  BYQ1BI2),  9V029I2I,  BYOUZI,  OTMIII, 
S0RV9(21 ,  9IV9I2I,  ORTA  (II.  RITAIII,  OrOlRIIl,  8T0IRIII 
DIMENSION  EL ( II. II I , EMI  12,11,  X111O0I,  M 1100 1.  13(1001,  14(1001, 
USI100I,  KOI  1001  .  XT  1 1001 ,  X9I100I X*(  100  i ,  X  10(100 1,  111(110), 
2112(100) 

COMMON  KCONI ,01, KC0N2,02 

EQUIVALENCE  (KCONI  111, LCONII  ,  (KCCHKIOl  I,  MCONI  I  ,  (KCON1IAOII, 

INC  ONI  I ,  I  KCONI  (901 1 ,  ICON  1 1 ,  (KCONI  11201 1,  ICON!  I,  (KCONtdl, 

2LC0N2I ,  IKC0N21 90 II  ,NC0N2I ,  IKC0N2(19Qll,NC0N2t,  (KC0NK2T01I, 
3IC0N2I  ,  (KCONI  ( 1401 1  •  JCONI)  ,  101(11, All.  (01(9011,911,  (01(9011, 

4C1I,  (01(4011,91).  (01 (1201 1, Ml  I,  102(11,92),  (02(401 1,021 , 

91031 19011 «C2I •  I02I2T01 I ,921 ,  (02(39011. H2I 

EOUI  VALENCE  (C(  11  ,E!I  ,  (EI2I  ,E2)  ,  (■(  II  ,E!>  ,( E(4I  ,E41 .1  EISl.tS)  , 

1(E  (91, 191,  (EITI,  «?l,  11(91,  E9),(B(4),E4I,(E(10),  1101,  (((III,  Ell) 
EOUI VAL ENCt  (NRO 1 , NROI , ( NTHl ,NtMI , INTI ,MT I 
■OUI VALENCE  (',011 , I J', KCONI I 

EOUI  VALENCE  (RHOtTITLEI  ,  ITHETA,  TITLE  191 1 ,  IMAM,  TITlEI  141 1 
C  INRUT 
2002  '0RMATIA9, 20131 

1002  'ORMAT  1 141 31 

1003  RORMAT  (114.7) 

1001  READ  INRUT  TARE  T.2002.RUN1 ,KOOEN,KTYR«,ROOSR,KCO  ,NAIOUT,KARR, 
lMSTAR,KARN,NtTAN.NACRCV,NRO,NTM,NMAJ 
NACTl-NACRCT«t 
TOL  -  10.44  (-NACYl; 

ILK  >  10 10 ( TOL  ,0 1 
tR  IKOOCNI  3494,4000,4001 
9444  RUN2  ■  4M  COM*  I 
RUN I  ■  4MNI0  NA 
0  RON*  •  492940909040 

00  TO  4002 

4000  9 UNO  >  9N  OILAT 
RUN9  ■*  9MATI0NA 
RUN*  •  AHL  NAVE 
00  TO  *002 

4001  RUN2  ■  AM  SMEAR 
RUNS  •  4M  NAVE 

9  RUN*  -  404040909040 

4002  MRIT6  OUTAUT  TARE  A. 4001, RUN1,RUN2,AUN3,AUN4,K00<N,KTYEE, RODS'. KCO 
1 . MAXOU T  •  K9"  , MIT  AR ,  KA RN ,  NS  T  AN, N4C9CT  . NRO ,  NT M, NNA J 

4003  RORMAT  1 12MIRUN  N(JN8«R  A9/// A9 ,A9, A4///T/4H  INRUT//10BM  KOOE 

IN  KTVRI  XOOSR  KCO  MAXTUT  RAM 

2  MSTAR  RARN  MST9N/9I 12//30M  NACRCT  NBO 

|  NTH  NMAJ  /4I12//I 

RE  AO  INRUT  TARE  T.IOOl.  <  II 1 1 ,1-1 .11) ,  1RH0I1 1,  I-l.NRO) ,  I  THETA!  I ). 
11*1  ?NTH)  ,  IXMAJI  D.lRl.HNAJI 
•4,’TE  OUT'*"  TREE  *, 1006, (f (ll.I-l.lll 


1008  FORMAT  (//5X.103H  H  /  R  B  /  COMT  NU 

1  NU  M  F  M  /  E  GAN  /  GAM  H/IP6E 18. 7//5 

2X»  85H8  ONEGA  /  CON  0M1  /  DM2  TAU1  /  ONI  TAU2 

3  /  0N2  C  N  TO  /  B/1P5E18.7) 

E7  -  E7  /  E2 

WRITE  OUTPUT  TAPE  6,t009,(RH0(  I  I  ,1-ltNROI 

1009  FORMAT  (///AH  RM0.5X.  IP3EI 8.  71 

WPITE  OUTPUT  TAPE  6,1010,( THETAf I ) ,I-lyNTM) 

1010  FORMAT  (//6H  TMF TA ,5X, 1P5E 18. 7/ 1 IX , 1P5E18. 7 ) 

WRITE  OUTPUT  TAPE  6,1011  ,( XNAJU)  ,  I«l  ,NNAJ » 

1011  FORMAT  I//11H  MAJOR  XBtR,5X , 1P5E18* 7/16X, IP5E18.7 > 

NMAJ1  ■  NMAJ  -  1 

>F  (NNAJl)  3001,505,499 

499  REAO  INPUT  TAPE  7,1002,1 INTERX( ( ) » I  *1 , NNAJl I . 

WRITE  OUTPUT  TAPE  6,1012 ,( INTERXt I i , I  •  1, NNAJl) 

1012  FORMAT  ( //7H  INTERX.5X, 10151 

500  XBAR(l)  •  XMAJI1I 
JAY  •  l 

DO  504  I  -1, NNAJl 

IF  ( INTERX t I ) )  3001,503*501 

501  INTER  •  INTERXd  ) 

ENX  «  INTER  *  1 

DEL  -  (XKAJU  +  l)  -  XMAJ(II)  /  ENX 
DO  502  J  •  1, INTER 
JAY  -  JAY  ♦  1 
YAY  ■  J 

502  XBAR(JAY)  •  XNAJU)  ♦  YAY*DEL 

503  JAY  •  JAY  ♦  1 

504  XBAR(JAV)  «  XMAJ(IM) 

505  NX  •  JAY 
FACT2  »  62 
EME*E(S) 

E(5I-ENE  *(ll.0tE(3ll/ll.0+E(4ll ) 

E6-EME*E6*( < 1.0-E4)/< I 1.0tE4)*IUO-2.0*E4l I )•( ( 1.0*631 
l*(l.0-1.0*E3)/(l. 0-631) 

E(6)*SQRTF(E6) 

62  -  F2  *  E6 

E( 1IM2.0-E (!))/( 2* OtElll) 

BETA-El 

BETAS0«E1*E1 

VMDUM-U.0-E4I/I  1.0-2*0*64) 

SVMOUM-SQRTF (2. 0*VHDUM I 
TR1  -  -1. 

TR2  •  0. 

2ETA  -  BETA 
KEY  *  -1 
JUMP  «  0 

IF  IKODENI  507,511,508 

507  KODEN  ■  0 
KEY  «  l 
REWIND  3 
GO  TO  530 

508  IF  (KEY)  509,1001,510 

509  KEY  ■  0 
GO  TO  530 

510  KODEN  •  1 
GO  TC  513 

511  IF  (KEY)  912,1001,1001 

512  KEY  «  0 
GO  TO  530 

513  REAO  INPUT  TAPE  7,10Q2,IEADX,KEYX, (KEYTHi I ) ,  !»1*NTH) 
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REAO  INPUT  TAPE  7, 1002 ,KTYPE ,KCO 
REAO  INPUT  TAPE  7, 1003, PHI ,AMPFAC,Ell 

WS  i  TE  OUTPUT  Tire  b  ,600?  , LCAua,KC YK  1 1  RE  YTHI  I  1 1  I  «  I, NTH) 

♦005  FORMAT  (35H1AODI  TIONAL  I  NPUT  FOR  COMBINED  WAVE//////21H  LEADX 
l  KFYX/I6, 115////6H  KEVTH/IX,  10 1 S > 

wSITc  OUTPUT  TAPt  6 ,6UU6 «KT YPE ,KCO 
♦006  FORMAT  I 22H  KTYPE  KCO/tlO,!12l 

WRITE  OUTPUT  TAPE  6 ,6007 ,PMI .AMPFAC.EI l 
♦007  FORMAT  I 56H  PHI  AMPFAC  C  M  TO  /  8/ 

I  IP  3f  18,  7  ) 

INTER  ■  LEADX  -  1 
IF  (INTER)  3001, 516, 914 

514  00  515  1  «  l, INTER 

515  LEADX  «  LEADX  *  INTERXi I  I 

516  XSTAR  -  XMAJ(KEYX) 

NT  -  NX  -  LEAOX  ♦  1 
INTER  ■  NTH 

J  ■»  0 

00  518  I  j  1 ,NTH 
IF  (KEYTH(I)I  517,518,517 

517  THETA! INTER*1 )  ■  -THETA! I ) 

J  *  J  ♦  1 

INTER  •  INTER  ♦  1 
KEYTH! J)  «  I 

518  CONTINUE 
NTHOLD  «  NTH 
NTH  »  INTER 

NEWTH  »  NTH  -  NTHOLO 
00  515  I  •  1, NTH 
.519  THETA! II  «  THE TA ( 1 1  -  PHI 
GO  TO  531 

530  NT  «  NX 
LEAOX  •  l 

531  A0J1  •  FACT2 
COUNT  -0. 

TALLY  -0,  _  -  . 

TAP  ■  0. 

IF  IKOOEN)  3001,533,532 

532  FACT2  »  FACT2  ♦  SVMOUM  ' 

TR 1  »  0» 

TR2  «  1. 

ZETA  «  “BETA 

533  TOT  »  E11AFACT2 
ADJ2  *  TOT 

ANGLE  *  TOT  -  2. *FACT2 
Ell  «  U  -  2*  *  TO  T 
IF.(XOOEN)  3001,537,536 
536  IF  (KEY)  3001,535,536 

535  ADJ1  *  AOJ1ASVMOUM 
GO  TO  537 

536  A0J2  •  A0J2  -  A0J1*XSTAR 

537  00  560  I  -  1, NTH 

560  THETA! II  •  .0176532927  •  THETA! I ) 

GO  TO  (8002,550,5501 ,KTYPE 
550  READ  INPUT  TAPE  7,1003,CTBIN,CAYIN 
EMIN  ■  CAYIN 
\  T1T0T-CTBIN6FACT2 
v  71T-T1*OT*TOT 
T2T  -  1.  -  TOT 
ACAPT-EMIN/d,  O-TITI 
GO  TO  (8002,80.10,8016)  ,KTYPF 
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8002  WRITE  OUYPUT  TAPF  6,8003 

8003  FORMAT  4///17H  RECTANGULAR  WAVE! 

<j<J  IU  t»v*u 

8010  WRITE  OUTPUT  TAPE  6,8011 ,CTB IN, CAY IN 

8011  FORMAT  4///16H  TRIANGULAR  MAVE///18H  C“P< TI-T01/8.17X, 1HK/IP2 

IFIB.71. 

GO  TO  8020 

8016  WRITE  OUTPUT  TAPE  6,8017 ,CT8!N, ENIN 

8017  FORMAT  4///23H  EXPONENTIAL  DECAY  MAVE///6X, l 2HC-M4 Ti-TOI/S, ITX, IHN 
1/1P2E18.7) 

C  CONSTANTS  FOR  L-NATRlX 
8020  PI  0.1615927 

PI  SQ*9. 8696067 
OiOUM»PI*Ei 

020UM«Q10UM*SORrF <42.0-2.0*^41/4  ».0“2.0*E3)» 

EL0YU-<2.0*F3) /( 1.0-2. 0«E3I 
CDUMM0.S-F6I/41.0-E6) 

008AR«E7**2*E8*<E9-I .33 33333 *CDUM*E10) 

01  8AP»E7*(  E9*EB*1.0“  li 3333333*C0UH*( E10*E8 )  I 
Q28AR* 1. 0-1. 3333  333*C0UM 
VOUM-4 1.0-F3l/(l,0-2.0*E3) 

YOUM1-OOBAR/PISQ 

YOUM2-E7**2/P1SO 

YDUM3-YDUN2*E8 

YDUM4»Q1 0AR*f 7*( E8+1.01/PIS0 

YOUN5-YDUM3*EB 

Y0UM6-Y0UM2«E10 

YDUM7«01HAR/PI 

YDUM8«E7*E  8/Pl 

Y0UM9.V0UM/VM0UM 

VDUM10»F7*< 1.0-f 401 /PI 

0360UM»PI*E2«SORTF(VOUM1/E6 

Y0UNH-E5/Y0UM9 

DO  9001  (ZR-l,NROl 

or  9ooi  izth«i»ntmi 

00  V001  UT-I.NT1 

A2 ( I ZR  • I2TH, I  £ T3 >0, 0 
LCON24 I2R, I2TS.I2T»»0 
82<I2R»12TH,I2T1 »0»  0 
MC0N2< I2R,I2TH,I2T*«0 
G24 I2P , IZTW>IZT)«0.0 
IC0N2(  !2R,I2TM,UT»«0 
H24 12R  *  I2TII,  I  ZTt  *0.  0 
JC0N2< !2R,I2TH,!ZTl-0 
C2(  :Zft,I2THTIZT)-0.0 

9001  NC0N24I2R,lZTH,IZT)-0 
C  MAJOR  100P  ON  P 

00  219  IP«1,KAPP 

GO  70  (9002,9003,9003) ,KTVPE 

9002  IFUP  2*UP/2>)  219,219,9003 

9003  P»IP 
PPI»P*PI 

COSST  •  COSF (PPI * ANGLE) 

SINST  -  SIUF (PFI * ANGLE ) 

JX  -  LEAOX 

00  l  I  TT»1 ,NT 

ANGL  -  A0J1*XBAR< JX)  *  A0J2 

JX  •  JX  »  l 

SINPT4 ITT)»SINF(PPI*ANGL) 

1  CnSPTUTT|-C05F<PPI*ANGL) 

C  CALCULATED  QUANTITIES 
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Ql«P*OLDUM 

01SQ*Gl*Ql 

Q2-P*Q/>DUP 

men  ^ f* r\t 

HALFSQ- ~^5  •  02  SQ 

PSQ.pop 


Vi  JO  t  *1  d  VOllM  t  /(If  A  «  UA1*  AA  i  *»  <  WIMIM«  #n«>  A«  «  a  «  .  « 

«*.  4  p  ■*  4  M  tfwnj/  rjKi-uvrv  f  IMAAf 


YVPJ-YVP1M  <YDUM5/PSQI«(  <  YD0*Z/PSQ)»1.0IHYDlJP&'P$Q|+1.0t 
YVPW.0*¥0UR*YVPl 


YSP 1*(  I  YDliWfc/PSU*  1«  01  /I YQUN2/PSQ+l*0l  J  *V0UM9 
VVP2* J Y0UA7/P1*(  Y0JMJ/PS0-l.0l-( Y0t*U/PSQ-Q2BARI 
1*<  Y0UW8/PCE7/PPI  I 

YVP2»YVP2/(  (  YDUR5/PS0I  *1  YDU*2mO*UO»  *YOUM2/PSQ*l*OI 
V VP2*2. 0*VOUApYVP2 

VSP 2*<  ( Y0UN10/P)  /(YOUA2/PSQU.O»  I*Y0UA9 
VP2-YVP2M. 33333333*YSP2 
VPl*YVPl4l. 33333333*YSPl 
SP-YP2**2/VPl**2 


SP*SQRTF ( 1, 0*SP| 

5SP*VSP2**2/YSP1 **2 
SSP-SQRTF<l.O*SSP> 

Q3BAR«$QRTF(  U*0«- SP>  /  YP  l ) 

Q3BAR«Q3BAR*P*Q340UN/SP 

03TIL*( i Y  P2/Y  P 11  /  ( SP+La 0 1  I *Q3BAR 

Q4BAK-S0RTF  (  <  1.0+SSP  I / YSPl  I 

04BAft»Q4BAR*P*Q340UA/SSP 


Q4YIL-1I  YSP^/YSPl  J/ISSPM.OI  >*Q4fl4R 
YVP 1 3U»YVP l*YOUAt 1 


YYP23U«VVP2«VOUN|| 
YSP12U-VSP1*Y0UH11 
Y$P22U*VSP2*YOU*ll 
ELRU«VVPi3l>-Q. 66 6666666* YSP12U 
ELTU-YVP230-0,  66  6666666+ YSP221! 
U8U-YSP12U 
UTU-YSPP2U 
EPDUPl “EL8U/ YOUP l 1 
EP0UN2*ELTU/Y0UP 1 I 
EPSB-FPOUHl*  YP1+EPOUH2  *YP2 
EPSB*EPSB/< YPl**2*VP2**2» 
S3*S0RTF  (Q3BAR  **2*Q3TIL**2 ) 

Si “SOP rFI04BAR**2*04TIL**2l 
03F-S3 


Q4E-SA 

ANGQ3*Q3TI L/Q38AR 
ANGQ4<*Q«TlL/04BAR 
PHI  3*- 4TANF(ANGQ31 
PHJ4*  -ATANf ( ANGQ4) 

If(KOOEN)  4 03 1*4031*4032 

4031  ABO-SQRTFIQ3BAR**2*Q3TIL**2) 

CALL  BESSEL  (K APN*A80*PHl 3*REJ*SMJ*PNGI 
IF(ENG)  4033*4035*3001 

4032  ABQ*SQATF(Q4BAR**2t94TIL**2t 

CALL  BESSEL  (KAPN*ABQ,PHI4*REJ *ENJ.ENG I 
IF(ENG)  4035  *4035*  3001 

4033  Q129-Q1*Q1*EI 
9182-Q12B*EI 
022B*Q2*Q2*6 1 
02B2riQ22B*E 1 

0320*0 3BARV42- 33 TIL**2 
04 20*04BAR**2-04T IL**2 


BSUH*0*3«EL8U40BU 


TSUW^U.  “.*tL  IU*U«  U 
Q3P“03BAR*Q3TI L 
Q4  p «Q4 BAP *04 T  I  U 

BBSUM«0.  ?*03F*IQiBAft*UBU4Q3TIL*UTU I 
BTD  IF^O,  5*Q3E*  <Q3BAR*UTU-Q3TI L*UBU) 

EPST-EPDUM2*  YP1— EP0UM1  *YP2 
EPST-FPST/I YPi**24YP2**2) 

5  00  6  IZOUT-l ,KAPN 

bjoib:  izoun*o.o 

BJ02B ( I ZOUT ) «0»  0 
8JQliIZOUT)-0.0 
BJ02I IZ0UT)*0.  0. 

BP  J  3(  I  ZOUT )  «C.>  0 
B’ J3( I ZOUT) *0. 0 
BP  J4(  I  ZOUT)  *0»  0 

6  6!  JMI  ZOUT)  »0.  0 

C  TILDE  FUNCTIONS  INDEPENDFNT  OF  ROt 

FNG-O. 0 

7  CALL  TILOE  I KAPN, "l .BETA, 0,0, 0JQ1B .HASTE , 

IBYQ1B. HASTE, ENG) 

TF(ENG)  8,8,3001 

8  CALL  TILDE  ( KAPN ,02 , BET  A ,0. 0 ,6 JQ2B , HASTE , 
lBYOrS.HASTF ,FNGI 

H22  «  ;.  *  BY02RI2) 

IMFNG)  9,  9,  3001 

9  C«!L  TILDE  (KAPN,Ql,l,0,0,O,0J0l ,HA5TF, 

16 Ywl, WASTE, ENG) 

IF(FNG)  10,10,3001 

10  CALL  TILDE  ( HAPN  ,02 , l. 0, 0, 0, B'Cr , WASTE , 

18Y02, WASTE ,ENG) 

IF  I  ENG  1  11,11,3001 

11  CALL  TILOE  (KAPN,S3,1.0,PHI3,BRJ3,BIJ3, 

1  Bit  Y  3  ,  B  I Y  3 ,  E  NG  ) 

IFIFNG)  12,12,3001 

12  CALL  TILOE  UAPN,S4,1.  0, PHI4 ,«R J4 , BI J4 , 

IBP  Y4.BI Y4.ENG) 

IFIENG)  13,13,3001 

13  NSTAR-0 
KPOE-1 

C  DO  LOOP  ON  ROE 

00  210  I«OE-l,NROl 

C  TIlOE  FUNCTIONS  OE PE NO I NG  ON  RUE 

RCE  -  RHOIIROEi 
1 A  DO  17  l  ZOl'T-l  ,KAPN 

BJUIRI  IZ0UT)>*0.0 

17  BJ«2R( 120UT!*0.D 

18  CALL  TILOE  {KAPN ,01 , ROE , 6, 0, BJwl R, WASTE , 

16Y0IR  .WASTE, ENG) 

IF(fNG)  19,19,3001 

19  ChLL  TILOE  I KAPN , 02 ,  ROE , 0. 0 , BJ02R , /AST F  , 

1BY02R  .HASTE, ENG) 

IF  I  ENG )  20,20,3001 

C  ZERO  CASE  FOR  ALL  THETA  ANO  T/T 

20  GO  TOiZl ,90311 ,KROE 

21  00  22  IEL-1,6 

DO  22  JEL-1,4 

22  EL(IEL,JEL)-»G.O 
IrlKODEN)  9022,9022,7023 

9022  EL  <1,  t  )*‘01SQ*<  (0.5 /El  I  *B  JOIB  <  2  »- <0.5*EL0VUM  .0  I  *B  JQ1BI  l )  ) 

EL  (l,2)-0lSQ*H2.O/012B)  *BYQIBI  ?  )  -  (0.5*EL0VUM.  0)  *BYQl  B<  1  )  ) 
FL  (2,  1)«OISQ*<0. 5*BJ01 (2 )- 10.5*EI  OVU* 1.0 l*BJ01 1  1 ) ) 
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EL(2,2l-01SCv(  <2.0/«lSQ»*BYQl<2>-<0.5*ELOVim.0)*aYQl<  1)  » 
EL231*BSUM*(Q32D4BRJ3I lt+2,0*Q3P*BI J3<  II  I 
k-7Sun.  iw 320* o i  ji  \  i  »-<•  jh  l »  j 

EL2  3.»«EL23i-  ( BB  SUM*BR  J3  <  2  I  -  B1 0 1 F  40  U3  <  21  I 
FL233-B5;tM*(Q320*B£  Y3I1I-2.04Q3P4BRV1I  II  » 
i*T Sir*!*  Y4(  n*2,  a*03P«>BI  Y3U) ) 

cl234-<2-0/33f  )* 1 Q3BAR* 1 UTU*8R Y3 12  I  ♦ll.'U*BIV3  (2)1  >Q3T  IL  * 
1(UTU*BIY3(2»-UBI'PBRY3I2>  I  I 
£L<2,3>*eL232+FL233-PL234 
Ft (3,1  I »(-0l SO  I *0.5*8401(2  I 
EL13,?I*(-2.0)4BYQ112) 

EL43’-0.  !*03E*1Q3BAR*BRJ312) ♦03TIL*B! J3(2I I 
EL332-12.0/Q3E  >*<GI8AR*BI  Y3  U> -Q3T I  L*L«Y3  <2  I  I 
FU3,3)«EL33l*EL332 

EL261=BSUM*(032D*BlJ3(l)“2.0*QiP*BRJ3( 111 
UTSUM*  1Q32D4BR  J3  ( 1 1»  2-  0*Q3P*B1  J31 1 1  ) 
EL261«EL261-(BTDIE*BRJ3{2l«-BBSUM4Bi  J312I  I 
EL  262*BSUM*  (  0320*BR  Y3  1 1 1  ♦2. 0*S53P*B1  V3(  II  ) 
[♦TSUM#<2.0*Q3P*aRY3m-0320*  ?i  Y3C1I  I 
EL26J--EL262 

EL263-12.0/Q3E  )4  <  Q3BAR4  (  UBU*BR  Y3  <2 1-UTU4BI Y  3  12  I  I  «• 

1Q3TI L* (UBU+B 1 Y3 ( 2  I ♦OTU*BRY312l  1 1 
EL(2,6  I •EL2M»EI. 262»EL263 
EL ( 2 ,6 )““EL (2,6) 

EL361-(-0.«i)*Q3E*((33BAR*BIJ3(2l-Q3TlL*BBJ3(2)l 
EL  362« 1  2*  0/C3F I* ( 03BAR4BR Y3 12)4Q3TIL*BIY312>  ) 

EL  1 3,6  I  «EL361*EL362 
FL15,3>~EL12,6I 
EL16,3)—EL1  3,6) 

00  23  IEL-1,3 

IIEL-IEL*3 
DO  23  JEL*1 ,3 

JJEL*JEL*3 

23  ELI l IEL* JJEL )“EL 1 IEL » JEL) 

CALL  2KHSP  (P,E2,E5,E6,Ell,03TIL,03BAR,EPSB,EPST,VMOUM 
l,KTYPF,KAPPfKOOSP,CAY?N,TlTOr,TOT,TlTfT2T,RE J. EMJ,  EH,  AUPT ) 
GO  TO  9024 

7023  FL  1 1  » l  )*(0. 5*02SQI *1 B J02B ( 1 >- 1 1.0/E  1 1 4BJC28 1 21 ) 

FL(  l,2l°0.5*Q2S346YQ2B< l)-(2.0/£L!*BYQ2BI2) 

EL ( 2 , 1 1 a ( 0» 3*02 SO  I *1 BJQ2 (II-  BJQ2 (2) I 
EL ( 2,2 l«( 0. 5*02SQ;p(8Y02  Ul- (4.0/0250) *BY02 121 ) 
EL231-Q4E*(04BAT*BRJ412 l+04TIL*BI J4<2) l-U420*BRJ4(l»'2.0P04P 
1 *B 1 J4(  ll-0420*Bl  Y4U  ) *2. 0404P4BRY4 (1) ♦ U. 0/04E  > 
2w<048AR*8JY4<2)-04TJL>?BRY4<2)  ) 

EL231-0.5*EL231 

El  2.s2»Q4f  4  ( Q48AR48 I J4 ( 2 )-04T IL*8R J4 1 2 )  I -04 20*8  I J4 ( 1 ) ♦2.0*04P 
1*SPJ4( 1 1 ♦Q42D*BR  Y41 1 1+2. 0*Q4P*Bl Y4 I l )- (4.0/04E I 
2* 1 Q4BAR*BR Y4  <2 )♦ Q4TI L*BI Y4 ( 2  I ) 

HL232-0. 54EL232 

EL  1 2 ,3  I «UBU*EL 23 l- JTU*EL232 

EL ( 3, 1  I  •  HALF  S0*8 JQ2 ( 2 ) 

EL  <  3,2  I  ■  2.4 BYQ212) 

£113.31*12. 0/04E )*( Q4TI L*BR V412 ) -04 BAR *BI Y41 2) > 

1-1  0.  9*Q*E  )  *1  Q4BAR4BR  J4  (2  )  4<34T  IL.4BI  J4  (2  I ) 
EL(2,6)-(-1.0)*1,jTU*EL2  314UBU*EL232» 

EL  36 l" 1 2. 0/04E )• <  Q48AR4BRY4 12 1 4Q4T l  L*B I Y4( 2 )  I 
1—1  0*  4*04E )  * ( Q*8Atl*BI  J4  12  I-04T I L*BR J4 (2  >  I 
EL ( 3,6 1 *-EL  361 
EL(9,3»  —  EL12,6) 

EL  <6,3I«-EL( 3,61 
00  7024  IEL-1,3 


IIFL*1EL>3 
00  7024  JFL-1.3 

JMri.- JcLtj 

7024  ELI  IIEL,  JJELI-EH  1EL,JELI 

CALL  ZRHSM  ( P, E2 ,E5 ,E6 .E 11 ,Q4T IL ♦ Q4BAR .EPS 8, ERST , VMDUM 
i  tKTVrS  ,KAPP  ,  KGOSF.CAVt  Fi,  TiTOT  ,  1  Ot  *  T  i  T  *  T  2  I  »KtJ»tHJ,tH*  ALAK1  » 

9024  IF  (  KCO  I  27,27,25 

25  WRITE  OUTPUT  TAPF  6 ,26 ,  (  (EL ( I , J I  ,  J» l  ,6  I, 1-1,6  I , ( EM( K ,  1  ) ,K»l , 61 

26  FORMAT  (36H1THF  MATRICES  FOR  N-0  ARE  AS  FOLLOWS  //// 

1I3H  THE  L  MATRIX  //6E17.7/  6E17.7/  6E17.7/  6E17.7/  6E17.7/ 

2(r 17.7///13H  THE  M  VECTOR  //  6E17.7/// » 

27  CALL  LINS7S  (EL, FM,6 ,LFROR I 
IF (LFRUR I  3006,3006,9027 

9027  IF(KCO)  30,30,28 

28  WRITE  OUTPUT  TAPF  6 ,29 , ( EM< K , 1 ) ,K* 1 ,6! 

29  FORMAT  < I3H  THE  X  VECTOR  //  6817*7///////! 

30  IF(KOOEN»  7031,7031,7030 

7030  XU  1)«0„0 
X2(  1!  -  0. 

X3(ll  *  EM( l , 1 ) 

X4< l)*EM(2f II 
X5(  i!«0.0 

XF(  1 1 “EM ( 3 , 1 1 
X7( 11*0*0 
XF(li  -  0. 

X9( 1)  -  EM(4,l I 
XIOC  l!  *EM(3 , 1’1 
111(11*0.0 
X12(  II *FH(6, It 
GO  70  7032 

7031  XI  (  1  > *FM (1,1! 

X2( 1 1 *EM (2,11 
X3(  11*0.0 
X4<  1J*0.  1) 

X5UI*EM<3,1) 

■XM1|.»0,0 
XM1  )*EM  (4. 1  > 

XP( 1  !*EM (5,1! 

X9 ( 1 ) *0*  0 
X10(l)*0.0 
Xll(U*|M(6tn 
X12( 11*0.0 

7032  KROF-2 

9031  IF(KCO)  9034,9034,9032 

9032  WRITE  OUTPUT  TAPE  4, 9033, ROE 

9033  FORMAT  ( 5H1R0E-E 14* 7 , 5X, 39HMATRI X-X3AR , SRR , STT ,SRT ,UR ,UT  FOP  N*0 
1//! 

9034  A  11*01  SO* ( 0*  5* ELOVU* l. 0 1 
A  12* 1* O/ROF 

Ali*Ail*(BJQlRU!*Xl  (1  UBYOIRI  1  I  *X2(  1)1-01  SQ*A12*(  BJO 1  R(  2 1  *Xl  (  1 1 
1*0. 5*< 2,O/OlSai*HV01R(2)*X2( II ) 

A14«All*(BJQlR(l)*X7(l)*BVQlR(lI*Xtt(ll »-QlSQ*A 12*( 8101 R ( 2 !* X7 ( l I 
1*0. 5* ( 2*  0/01S0 1 *8 V01 R ( 2 ! *X8 ( 1! ) 

Bll  •  .5  *  ELOVU  *  01  SO  *  ( BJ01 R ( 1 ! *X1  ( 1 )  ♦  BYQIR( l»*X2( 11) 
B12*A12* ( (Q 1SQ/2.0 )*BJ01R( 2) *X1 1 1 1 *2.0 *8 VOIR (2  J *X2 ( 1 > I 
B13  -  .5  •  ELOVU  *  QISO  •  (BJQ1RI  ll*xmi  ♦  BVQ1R  ( 1  !*X8(  l ! ) 

B 14*A1 2* ( (0  ISO/2. 01 *6 J01R (21 *X7 ( 1 )t2.0*8VQlR(2)*X9( 1 )) 

Cli  ■  8J02P ( I )  -  8J02R (2  I /ROE 
C 12  »  8Y02R( 1)  -  4,*eVQ2R(2»/(02Sa*R0F » 

Cl 3  •  HALFSQ  *  (C  1 1 * X 3 ( 1 )  ♦  C12*X4(1») 

C 14  -  MALFSO  *  ( C ll*X9( 1 »  ♦  C 12* X 10 (  II  ! 
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Gl  la(QlSQ/2. 0)*UJQ1R(2)*XI (1 )42.0*BYQ1R< 2>*X2(  1) 
Gl2»(QlSQ/2.  0) *B  JQ1R (2  >  *  XT  1 1 ) +2. 0*BYQIR ( 2 ) *X8( 1 ) 

ni*  -  nilr  av  ♦  o  J  w2  M£l 

H12  «  *  BYQ2R  ( 2t 

H13  -  HU  •  X3  ( 1 1  4  H12  *  X4  (1) 

H14  »  r>  i  1  *  X9(i»  4  HI2  *  Mi  0 1  i  i 

H21  ■  HALF  SO  *  63028(21 
H23  ■  H21  *  X3( l I  ♦  H2 2  *  X4(ll 

H24  -  H21  «  X9(l»  4  H22  *  XlOtll 

H25  »  H23  *  CQSST  -  H24  *  SINST 

00  44  I TH  »1 1 NTH1 

33  IF  (KCOI  36,36,34 

34  THFT  «  57.2957796*THETA(ITH( 

WRITF  OUTPUT  TAPE  6,35,THET 

35  FORMAT  <////  7H  THETA-E14.7  tf\ 

36  KSTART-1 
KSTOP-6 
KWASTE-t 
IND1  *  ITH 

IND2  «  3* I TH  -  3  4  IROE 
DO  64  I TT  *1 ,NTi 

IF(NYl-ITT)  37,37,38 

37  KST0P«NTl 

30  *ASTE(KWASTE I  •  XBAR(ITT) 

K W AS TF  aK WASTE* 1 
DO  410  INK  •  1,5 
KON  •  KC0N2IIN02) 

IF  (KDN)  401,402,402. 

401  D 1  (  I  NO  1 1  -  0. 

GO  TO  409 

402  GO  TO  <403,404,405,406,4071 , INK 

403  A15-A13*COSPT< l TT I -A  14* SI NPT ( I TT I 
OK  IND1I  -  A15  *  TR1 

GO  TO  408 

404  B15><B114B12I*C0$PT(!TT|-| B13*B14I*SINPT< ITT  I 
Dl(INDl)  -  615  *  TRl 

GO  TU  4CB 

405  0  5  ■  C13*C0SPT  (ITT)  -  C14*S  INPT  (ITTI 
01IIND1)  -  C15  *  TR2 

GO  TO  408 

406  G15«Gll*C0$PT< ITT  I -G 12* SI NPT ( ITTI 
OK  IND1I  •  G15  •  TRl 

GO  TO  408 

407  H15  *  H13  *  COSPT(ITT)  -  H14  *  SI NPT 1 1 TT )  -  H25 
Dill NO l I  ■  H 15  •  TR2 

408  KC0N2I IND2 I  -  1 8  I G(01 < iNDl  I  , KON I 

409  INOl  -  INOl  4  300 

410  IN02  ■  I  NO 2  *  900 

INOl  «  INDl  -  1490 
IN02  •  IND2  -  4470 
IF(KWASTc-7>  9039,40,40 

9039  IF  < KSTOP-I TT )  44,40,44 

40  KWASTF-1 

IF (KCOI  44,44,41 

41  WRITE  OUTPUT  TAPE  6,42, «WASTE( III ,11*1,61 

42  FORMAT  (5X.6fl7.7I 

WRITF  OUTPUT  TAPE  6,42, (All  ITH, 121 , 1 2* K ST  ART , KSTOPI 

WRITE  OUTPUT  TAPE  6,42,(811  ITH,I3),|3«KST ART ,K$TOP| 

WRITE  OUTPUT  TAPE  6,42,(C1(  ITH, 141 , I4-KSTART, KSTOPI 

WRITE  OUTPUT  TAPE  6,42,(G1(  ITH, 15 1 , I5-KST ART, KSTOPI 

WRITE  OUTPUT  TAPE  6,<*2,(H1(  ITH,  16  I ,  16*KST  ART, KSTOPI 
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K  STAR  T-K  STAR  T«-6 

KSTQP-KSTART+5 

HR I  TP  OUTPUT  TAPE  6,43 

43  FOPMAT  I /// » 

44  CONTINUE 

00  45  I ZER-l ,  NTH1 

DO  45  JZFR«l,NTl 
LC0N1 (  I ZFR,  JZER 1*0 

MC0N1 (  !ZER,JZER)-0 

NCONK  rZER,JZER»«0 

(CONK  IZER,JZER)»0 

45  JC0N1 t  I ZER, JZER ) -0 

NSTOP-KAPN 

E  33-1.0 

F36-1.0 

LOOP  ON  N 

00  140  N«l,NSTOP 

IFIKOOEN)  9047,9047,9046 

9046  FN«-N 

GO  TO  9048 

9047  FN«N 
904 S  EN-N 

Nl-Ntl 
N2-N42 
ENl-EN+l. 0 
ENM1-EN-L.0 

IF(N-NSTAR)  59,59,46 

46  BYD1-0YOIB( 2 ) 

6Y01B(2l-(EN/<ENt»BETA)l*BYQlB(2l-U0.25*01SQ)/(EN«ENl))*BYQ14m 

BYCIBI 1I-BY01 

BYD2-8VQ2BC2) 

RYQ2B1 21 -I EN/( ENI*BETAU  *BY02B(2 1“ ( ( 0. 25*Q2S Q) / ( EN*ENl ) I •BY02BI l ) 
BY02BI II-8YD2 
9030  BYD3-BY0LI2) 

BYQ1 ( 2 1 ■ ( EN/EN l) *BVQ1 (2I“ I  0«25f ( EN*ENl 1 1 *BYQl( 1 1 *Ql$Q 
BY01 ( 1 ) -BY03 
BY04«BYQ2<2» 

BVQ2 ( 2 ) * ( EN/ENl I  *BYQ2(2l-l0.25/(  EN*ENl  n*BY02«U*«2S0 
BY02 ( 1 1 -BY04 
BY05-BRY3I2) 

BRY3(2»-(03e/S3»*(EN/ENl)*(C0SF(PHl3)*BRY3(2)7SINF(PHI3)*BIY3(2) I 
l-(  <0.25*03E**2I/IEN*EN1H*BPY3(U 
BRY3UI-BY03 
BYOE-R I Y  3 1 2» 

8IY3(2)-(Q3F/S3>*IEN/ENIJ*{C0SF(PHI3)4BIY3(2I-SINP(PHI3)*BRY3(  II) 
L-H0.2  5*03E4*2>/(EN*EN1I  )*BIY3(1I 
8IY3m-BY06 
BY07-BRY4 I  2 1 

BR  Y4( 2 ) • <  Q4t /S4I * IEN/ENL I *ICOSF ( PHI4 I*BRY4(2 1+S INFIPHI4) *BI Y4(  21 ) 
l-(  (0.25*Q4E**2)/IEN*ENin*BR Y4UI 
BR  Y4( 1 ) -BY07 
BY0B-BIY4I2) 

BI Y4<2I<<04E /S-,1 *(EN/ENII*(COSF I PHI4I-BIY4C 2 l-S INF(PHI4)*RRY4( 1)» 
l-(  (0*  25*046**2  I  Z(FN*ENIH*BIV4<1» 

BI Y4(l»-BY06 
COMPUTE  THE  l-NATRIX 
ELMl.0.5*ELOVU7l.O-EN*FNMI/01B2 
EL M2- ( 2. 0*EN*E NM 1  /Q2B2  > -I. 0 
ELM3-0,5*EL0VU+1.0-EN*ENM1/01S0 
EL  M4»  ( 2.  0*E  N*E  NM1 7 02  SO ) - 1 »  0 
FJ5-F35*(03F/(2,  0*tNM 
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F36-F36*(Q4E/(2.0*FN)  1**2 
EL t  l.ll-ElMl*BJQia(Nn-SJ9lB| 


)/<2.0*El*ENll 


FU  |,2)-ELM1*BYQ1B<  1I-BYQIBI  2  ) *2.0*FN1/QL2B 
ELI  I *2  )■  !-QlSO l*EL ( 1 f 2 ) 

EL{  l*3l=EN«l*3J32Bt  Nii-ij22t»*#jy*D<N2l/<2.0*ENl» 
El<l,3)«<FN/BETAS0)*EL(l  *31 
EL(lf4)-ENHl*RY92B(  1I-2.0*E1*EN1*BYQ2B(  2) 
Fl<l,4»»<FN/BETASQ>*EL<l.4l 
00  4?  KM»1,2 

00  47  LM»5,12 

47  EL  (  KM»  LMI  »0.  0 

EU2,ll«Q12B*BJQlB<  N2  ) /(2.0*ENI  I-ENH1*BJ91  B(  Nl I 
El <2,1 l«(FN/BETASQ)*EL<2.1) 

El(2,2>-2.0*El*ENl*8VgiB<  2  I-ENMI*BYQ18<  1» 

FL(2,2»-(FN/BF  TASO>*EU2,2l 
E  l  (  2 « 3  )  ®E  H2*B  JQ2B(  Nl )  +B  JQ2B  (  N2  >/<El*ENI) 

£U2,3l«<-0.  51  *Q2S9*EL(2  .3) 

EK2,4I»EIM2*BYQ2B<  ll*4.0*ENl*BYQ28l  2  I/022B 
El ( 2, 4 1« < -0.51 *02 SQ*El 12, 4) 

£1(3,1 )»ELM3*BJ01 (Nl I-BJ01 (  N2  )*0.5/ENt 
61(3.11  —  Q1SQ*EL(3,1I 

EL(  3,2I«ELM3*BVQI<  ll-BYQM  2  I *2.0*ENl/QIS<i 
ELI3.2  I— QISQ*EL(3,2I 

El  4  3, 3  I »ENMl*B JQ2 <  Nl )-Q2  SQ*B JQ2 (  N2  l*0.S/ENl 
EL(3,3)-FN*El(3,3l 

El<3,4 l«FNMl*BYQ24  l 1-2, 0*EN1*B YQ2 (  21 
ELI3,4|«FN*FU3,4I 

FL35“B  SUM* (0320* BRJ3 (Nl I ♦2.0*Q3P*Bl J3( Nl I I 
L-TSUM*  <0320*61  J3(Nl»-2*0*Q3P*BRJ3(Nin 
EL35i»EL35-(BBSUM*8RJ3<  N2  I -BT0JF*B t J3(  N2  11/ EN i 

1- EN*ENM1*4UBU*BRJ3(NI)-UTU*BIJ3(N1) I 
EL351»F35*EL351 

FL352«8SUM*( 0320*81 Y3(  1 )-2. 0*03P*BR Y3 (  111 
l+TSUM* (Q32D4BP Y3 (  1 1 *2. 0*Q3P*Bl Y3 (  111 
E1353-EN*ENM1*(URU*8IY3<  ll+UTU*BRY3<  1 1 »♦< 2.0/Q3E * *FN1 
1*(Q3BAR*(UTU*BPY3(  2I*UBU*BIY3<  2 1 I+Q3T IL*(0TU*8IY3(  21 

2- UBU*BRY3<  2)11 
El(3,5l»£L351*EL352-EL353 

EL36l«ENMl*BRJ4(Nll-0.5*<04E/ENll*(04BAR*BRJ4<N2l 
1*Q4T IL*Bl J4( N2 )  I 
EL361»F3t*EL361 

EL362-ENMi«BIY4(  H-2.  0*  ( ENI /Q4F  I*  <  Q4BAP*BI  Y4(  2) 
1-Q4T!L*BRY4<  211 
RRRS-EL361+EL362 
RPRS*~RRR5 

EL363»ENMI*B1 J4< Nl I -0. 5*( 04E/ENI I  * (04BAR*B I J4< N2 I 
l-Q4TIl*BRJ4(N2>) 

EL363-F36*EL3E3 

FL364-ENMl*BBY4<  l 1-2. 0* (EN1/Q4E )*<Q4BAR*BRV4<  21 
I+04T I L*BI Y4(  2)1 
EIRRS-EL344-EL363 
EL(3.M-UBU*RRRS-UTU*E!RRS 
El <3,A)-FN*EL(3,6I 
00  46  KM-3,6 

:  00  48  LM-7,10 

48  FL(KM,lMI»0.0 

FL  95 »B SUM* ( Q 32 0*81 J3 ( Nl  1-2.0*Q3P*BR J3<  Nl  1 1 
l*TSUM*(032D*BRJ3(Nl)*2.0*Q3P*B( J3INII > 
EL95l-EL95-(8T0IF«dRJ3(  N2  )*BBSUM*BI J3(  N2  ) I /ENI 
1-EN*ENN1*(UBU*BI J3 ( Nil tUTU*BR J3( Nl 1 1 
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EL95l<F35*EL951 

EL952<BSUX<<O32D<0RY3<  1 )<2. 0<Q3P<BIY3<  111 
XTSUH* (2. O<O3P<0RY3<  1 1 -0320*01  Y3 1  11) 

Pi 

EL953<EN<ENH1<<UBU<6RY3<  1)-UTU<BIV3<  IK  <  <  2.C/Q3E )  *EN1 
1<<O30AR<<U0U<BRY3<  2 l-UTU*Bl Y3 <  2II<Q3TIL<<UBU<91Y3<  2) 
r*nrii*«Dv j i  2)1) 

~FL<9, 5<<EL951<EL  952<EL953 
ELi3,  U)«-EL<9,3) 

EL  <  9f ft  I <UTU<PRRS<UBC<E I RRS 
EL<9,6)-FN<EL(9,6< 

EL13.12I—  EL(9,ft| 

El (4, l<< <Q1S0<BJO1<N2 1 l/(2.0*ENl )-ENHl*0JOl< Nil 
EL  <4t 1 l*FN<EL (  A* 1 1 

EL<4,2i»2.0<ENl*BYQH  2)-ENHi<SYQl <  II 
EL<4,2I<FN<EL<4,2< 

EL<4,3X-!Q2SQ/2.  0)  <<ELN4*n JQ2< N1  > <QJQ2< N2  XENl ) 

EL  ( 4,4  X  ”  (Q2S0/2.  0<  •  (ELN4*8YQ2  <  11 ♦< 4.0/ 02S0 KEN L*BY02<  2)1 
EL451<ENMI<8RJ3< Nil-0. 5<<Q3E/fNl) <<Q3BAR<BRJ3<N2I 
1<03TIL<B<J3<N2 II 
EL451-F35*£L451 

EL452<2NM1<BIV3<  1 1 -2«  0*< EN1/03E 1<<Q38AR<B1Y3<  2) 
1-03T1L*BRV3<  2  1 1 
RR  TU<EL451<EL452 

EL433<ENM1<8U3<  Nl!-0.  5*  <Q3E/EN1)<  <O3BAR<0IJ3<N2< 
1”O3TIL<0RJ3< N2 11 
EL453-F35*FL433 

EL454<ENMI<BRY3<  l*-2.  0<<EN1/Q3E  X<  <J3BAR<8RY3<  21 
1<03TIL<BIY3<  211 
E1RT0-EL453-EL454 
EU4,5)-UBU<RRTD-UTU<EIRT0 
El<4,5)<FN*EL(4,5l 

EL461- <2.  0*£N*ENMl-Q420<*BRJ4<NU-2. 0<04P<BIJ4<NII 
1<<0AE/ENLI  M048AR<8RJ4<N2X04T1L<8IJ4<N2M 
EL*6l-F36*EL461 

EL4ft*<<2.0*EN*ENMl-Q420)*B!V4<  U<2<0<Q4P<8RY4<  11 
l<<  4.  0'*ENl/04E  I  *<  04  BAR  <B  I  V4  (  2XQ4TIL<BRY4<  2)) 

RRTS*  I<5<<EL461<EL462< 

EL 463< < 2. 0*EN*ENM1-Q42DI*B!J4<N1< <2. 0<Q4P<BRJ4<NI » 
K<O4E/ENl<*<O4BAR*0I J41N2I-04T <L<BRJ4<N2 K 
EL  46 3«F36<EL463 

EL464- <0420-2. 0*EN*ENMI  »<8RY4<  1X2<0<04P<BIY4<  II 
l- <  4. 0*ENl/04E<  <1 04BAR*BRY4(  2<<Q4TIL<B1Y4<  21) 

EtRTS*0. 5<<EL463<EL464< 

EL (4.X <U8U*RRTS-UTU*EIRTS 
EL  <10, 5<<FN<<UTU<RRT0<UBU<EIRT0 1 
EL (4, 1 1 1 »-EL  < 10* 5  I 
;L  < l 0, 6<<UTU<RRTS<UBU<E 1  RTS 
EL  <4, 12 ) ■— EL ( 1 0, 6 < 

FL(B, 1<<EN<BJ0KN<  <- t 01 SQ/EN1 ) *0. 5*B JQ1 ( N2 ) 

EL <5, 2<<EN<6V0H  1 )-2.0<ENl*BY0l <  2) 

ELI?, 3<<FN<BJC2<NI< 

FL<5, 4<<FN<8Y02<  U 

EL551<EN<0RJ3<N1 <-0. 5< < 03E/ENI ) * < 038AR*BR J3 <  N2 > 

1  <0  3TIL<6IJ3<N2K 
EL551-F35*EL551 

FL<52<EN<BI Y3(  l<- <2.0*EN1 /Q3E ) • < 03BAR<BI Y3 <  2> 
1-Q3TII.*BRY3<  211 
EL  <5, 5  <»<-!.  OK<EL351<EL552< 

EL  <  5,6  (<<-F N<<<F36<BRJ4<NI  XB1  Y4<  1 )  I 

6L5111 «E N< BI J3(Nl)-0. 5* ( 03 E/ENl!*( 03 BAP<B1J3<N2< 


l-Q3TIl*BRJ3<N2>> 

E  L  5 1 1 1  «F  35*E  L5  i  1 1 

EL51l2a(2«0*ENl/Q3E )*(03BAR  *BR Y3 (  2lfQ3T (I •Rtvli  »i» 

1-EN*BRY3<  1) 

ELI5,ll)aEL5lll«EL5112 
ELI5»12)aFN*(F36*BIJ4(Nl ) -BR Y4 (  111 
t-L(6«l)a(»FN|*BJQHNll 
EL(6,2)a(-FN>*BY0t(  1) 

El(6,3)a(-EN)*BJ02(N1)«-0.5*(Q2SQ/EN1I*BJQ2IN2) 

EL<6,4)aI-EN>*BYQ2I  11  >2. 0*ENl*BYe2(  2> 
FL(6,5JaFN*(F35*BRJ3INll«-B]|V3<  1)) 

FL6 6  l«EN*B« J4< NII-O, 5*< Q4E/ENI ) * ( Q4BAR*BRJ4{ N2 ) 
l*Q4TIL*BIJ4<N2>) 

EL  66 1«F 3 6* EL 66 1 

EL662aEN*BIY4<  1 )-2. 0* < EN1 /04E t * I Q46AR*BI Y4 <  2) 
l“04ftL*BRY*(  2)1 
EL(6»6laFL661*EL662 
EL  (6-  1 1 )  «FNMBRV3{  1  )-•=  35*91  J3<  N1 } ) 
FL6l2IaEN«BUA(Nl)-0«5*<04E/ENlH'<Q4BAR*8IJ4<N2t 
1-SHTIL*BPJ4(N2I> 

EL6121»l“F36)*EL612l 

EL6122«f N*BRY4 (  l)-2, 0* IEN1/Q4E >*« Q48AR6BRY4 (  2) 

1+C4T !L*B I Y4(  2)) 

EL ( 6i 12  I aEL612 14EL6122 
00  51  KH-7,12 

KM6«KM-6 
DO  49  LMal,6 

LM6lalN*6 

49  El <KM,LN)a-EL<KH6,LM61l 

00  50  LM-7,12 

LM62-LM-6 

50  El(KM,LM)aEl<KM6,lM62> 

51  CONTINUE 
NSTAR-N 

IFIKODENI  9051,9051,9052 

9051  CALL  RHSP  <P,E2,E5,E6,EIl  ,03TIl,03BAR,Ei>SR,ePST,VN0UM 
l,KTYPE,KAPP,KaOSP,CAYIN,TlTOT,TOt,TlT,T2T,REJ,ENJt ACAPT,N,E»M 

GO  TO  9053 

9052  CALL  RHGM  ( P,E2 , E5, E6 ,EU ,04T1 L.Q48AR ,EPSB, EPST , VMDUM 
l,KTYPE,KAPP,KODSP,CAYIN,UTOT,TOT,TtT,T2T,REJ,EMj,AC4PT,H,EM)  . 

9053  IF(KCO)  55,55,53 

53  WRITE  OUTPUT  TAPF  6 ,54 ,N, I(EL( I , J5 , J-l , 12  I , I -l , 12 ) , ( FN( K, 1 1 
1,K«1,12) 

54  FORMAT  ( 20H1  THE  MATRICES  FOR  N»I4,6H  ARE  ///// 

113M  THE  L  MATRIX  //6E17.7/  6E17.7//6E17.7/  6E17.7//6F17.7/ 

26E 17« 7//6E 17,7/  6E 17. 7//6E 17. 7/  6E17. 7//6F17.7/  6E  17.7//6* 17* 7/ 

36F 17*  7//6El7«7/  6E17.7/V6E  17.7/  6E17.7//6E17.7/  6E17.T//6EIT.?/ 
46E17.7//6E17.7/  6E17.7/////13H  THE  M  VECTOR  //  6E17.7/  6E17.7////) 

55  CALL  LINSYS  (FL, EM, l 2 .LERORI 

IFILEROR)  3006,3006,9055 

9055  IF(KCO)  58,58,56 

56  WRITE  OUTPUT  TAPE  6,57,  (EMU,  1)  ,K«  1 , 12  ) 

57  FORMAT  ( 13H  THE  X  VECTOR  //  6E17.7/6E17*?  ///////) 

58  00  9058  IMT-1,12 

IF (EMI IMT , l ) )  9059,9058,9059 

9058  CONTINUE 
NSTOR«N 
GO  TO  160 

9059  Xl(Ni)«EM( 1,1) 

K2(N1I*EMI2,1) 


X4(N1>-EMI4*1) 

XS(M‘*EM(5,  I) 

U  ,'Nl)*E*U,ll 

V<  (IMII  •*t«l  I  ,  i  ] 

X8HU)*EMI8,l» 

X9INn»EMI9„tl 

XI 01 Nl I  *FM( 10(11 
XI t(Nll«EM(ll, 1) 

XJ2(NU*ErtU2,l) 

IF  IN  -  U  904 «♦  9049, *59 
9049  Y01B1.  »  8 YQ1BI 1 ) 

YQ 10  2  -  BYQIBI2) 

VQ2B1  «  0YO2BID 
Y0282  «  BV028I21 

fl  ■  B  JQ1  B  (  2  )/8£  TA  -  01*0*83011(3)74. 

T2  »  YQ131/3ET4  -  4,*YO102 
T3  *  RJ02BI2 J/2FTA 
T4  -  Y02B1 VZET4 

TA  a  nniU)  *  T2*X2I2I  ♦  T3*X3 ( 2 f  ♦  T4*X4(2) 

TB  *  Tl*X7(2)  *  T2*X8(2)  ♦  T3*X9<21  <*  T4»X10<2> 

tG  *  TA*COSST  -  T6*S1NSV 

.  71  «  BJQ16I2I/2ETA 
T2  *  V01B1/ZETA 

.  T3  -  BjQ2B(2)/BETA  “  Q2SQ*BJQ2313>/4. 

V.  T.4  •  YQ261/BETA  -  4.«Y02,B2 

TA  »  Tl*XH2)  +  T2+X2I2I  ♦  T3*X3<2)  ♦  T4*X*12) 

TB  «  T1*X7<2»  ♦  T2*X9J2)  *  T3*X9<2>  ♦  T4*X10I2) 

■ ...  TH  ■  TA*CO$ST  »  T8*$IN$t  . 

59  IF  IKCOI  .  ■■-62,62,60 

60  WRITE  OUTPUT  TAPE  6, 61, POE, N 

.  61  FORMAT  <  5H  1R0E“F  I4«7 , SX, 37HNATR  i  X-XBAR  ,  SRR  »  STT  ,  5RT  ,UR  »(JY  FOP  N* 

lt*//l 

62  BY09-8YQIRI2I 

8Y0 tR ( 2) “ (EN/I  ENl*RQE  )  )*BYQIR  (.2  )"■  I  <0.25*Ql,SQ  >/ <  £N*£Nl  )>*B*2»U  1) 

OYCllP  ( 1 )  »BY09 

BYOIO«»V02P(2) 

BYQ2M2)«IEN/IicNl*ROE)>*C>ru2R<2)-<  <0.25*Q2SQI/I tN»ENl I  >*bVQ2R<  1) 
BV02RIII-BY0IO 

TA A«0»  5*E10VU- (E  N*ENM1 1/ IQlSQ*ROE*ROE )  *1.0 

TA l»Ol $Q*(  TAA*BJ  Q1R  I  Nl  )- (0.5 /  I  EN1*R(S£ )  ) *BJ01R (N2 !  I 

TA2-QISG*(7AA*BYQIR(  ll'-((2»Q*ENl  )/ (Ql$Q*ROE  1 1  *BYQlR{  2I> 

TA  3«(FN/RQE**2)*  ( ENMl  *BJ02R  INI I- (I  025Q*RQE  )/(2,Q*£Nl  I )  *BJQ2P  I N2  )  I 
TA4»(FN/R0E**2)*(ENMl*BVQ2R<  II-2.0*EN1*ROE*&Y02M  2)1 
VABorAl*  XI (Nl ) *T  A2*X2( Nl l-T A3*X3 INI )-T A4 *X4 1 Nl I 
TAOTAl*X7<Nn*TA2*X9(Nl)*TA3*X9(Nl)-TA4*X10tNl  I 
T0B»O.5*E(.OVU*QlSQ*EN*ENMl/RaE**2 
TBI»TB.B*BJCIIR(N1  >*(Q1SO/(2.0CEN1*ROE)  l*BJ(UMN2) 

TB2"TBB*BYQIRI  1  )*<2.0*ENI/»0E)*BYQ1M  2) 

TB3»(FN/fiaE**2)*(ENMl*BJQ2R<Nl)'I<Q2SQ*RQF)/(2.0*£Nl M*BJ02B(N2) I 
TB4*-(FN/RQE**2  I*  (ENMl  *BYC2Ril  1  l-2,0*ENl'-‘ftOE*BY02M  211 
TBC«TBI*XIIN1)*TB2*X2(N1 >*TB3*X3(Nl)*TB4*X4(Nl! 

TBl)»TB l*X7  < Nl )  ♦(“  82*X8  (Nl  )*T83*X9(NI)*T04*X10(N1  t 
01  SR  ■  QISQ  *  POE 

TC1«(FN/RQE**2)*((QISR/«2.0*EN1I ) *0JQ1R( N2 )-ENMl *B JOIR ( Nl ) ) 
TC2«<FN/ROE**r>*<2.0*ENl*ROE*8YO!RI  2)“ENMl*BY01R(  111 
TC  3*1  0,5  *02  SO  I*  I  <2,0*EN*ENMl/(02S3*ROE**2n*BJ02R<Nl  I-BJ02RIM) 
l>(  1.0/It  Nl*ROE  )  I  *8JslZ‘MN2)l 

TCA«(0.5*02SQ)*( (2.0*EN*ENM1/<G2SO*HOE**2I )*eV32R<  II-BYQ2RI  l) 
1*1 ( 4»  0*E  Nl I /( 02SQ*R0E 1 1 *8Y02R(  2)1 
TCA»TCl*Xl(Nl)*TC2*X2(Nl I-T* 3*X3 (Nl 1-TC4 *X4 < Nl ) 


TCB«TCl*X7(Nl)*TC2*X8<Nl)-TC3*X9INl)-TC**XI0<Nl) 
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9104 

9105 

9106 
910? 


65 

66 

67 

68 


:074 


TG1«(EN/RC/FI*8JQ1RINI  >-IQlSQ/(2.0*ENl)  )*BJQtR(N2l 
TG2*IEN/R0E I *BYQ1R ( 1 )-2. 0*EN1*BV91R( 2) 
TG3»<FN/ROE)*BJQ2RlNl 1 
YG4«(FN/RQE  )  *8Y02R  U) 

TGA-TG 16X1  (Nil  ♦TG24X2INl)6TG3*X3tNl)6TG4*K4«  fui  » 
TG8»TG1*X7IN1I  +TG2*XA( N1  (♦  TC4*if« f mi  UTutviiwm  > 

TH1.(FN/R0E)*9J01RIM  )  . 

TH2»(FN/RCE ) ♦BYQiR ( l ( 

TH3» (EN/ROE  I  *3  JQ2R  (Nl  )**  (  (02  SO  )  /  (  2«0*ENl )  i*BJQ2R(  N2  I 

TH4«( EN/ROE  I *8  YQ2R ( 1 ) «2»  0*EN1  •3TQ2RI2)  ’ 

TH6«THI*X1INI|6TH?*X2INI  |6TH3*K3<NI)6TH4*X4IN1> 
TH6«TH1*X7 (N1I*TH2*X8(NIItTH3*X9(NI)+TH4*X10(N1 ) 

DO  9107  ITT  *  l iNTl  ’ 

TAOIITTJ  •  TAC*S I NPT ( 1 TT I 


TflO*S INPTI I TT I 
TCA*C  OSPTI ITT ) 
TGA*COSPTI ITT) 
9104,9104,9105 
TGC(ITT)  -  TG 
TH5*C0SPT( ITT ) 
9106,9106,9107 
TH7IITT)  -  TH 


T8E ( ITT) 

TCC(ITT)  . 

TGC(ITT)  - 
IF  IN  -  II 
TGCIITT1  » 

TH7< ITT)  ■ 

IF  IN  -  II 
TH7IITTI  « 

CONTINUE 
00  129  I TH«1 , NTH! 

THET  »  THETA (174  I 
THAT  •  57, 2957796* THE  T 
IFUCOI  63.68,66 
WRITE  OUTPUT  TAPE  4, 67, THAT 
FORMAT  t////  7H  THETA-E14.7 
KSTART-1 
KSTCP-6 
KWASTE-1 
IF  I KOOEN ) 

ENDUM-N 


TAB*COSPT( ITT  > 
I 8C*C0SPT (ITT) 
TCB*SINPTI ITT) 
TGB*SINPT( ITT) 


-  TH6*SI NFTJ ITTI 


n\ 


9074,9074,9075 


TR I Gl "CO SF ( FNOUM * THET  ) 

TR IG2«$INF(EN0U*UTNET  I 
GO  TO  9076 
5075  ENOUM«-N 

73 IGl*S! NF I ENOUMYTHET  ) 

TR IG2»C0SFCEN0UM*THEr  I 
9076  IND1  *  I TM 

1ND2  «  3*1 TH  -  3  +  IR&E 
00  129  ITT»1,NT1 

C  CONVERGE NCE  TESTS 
00  110  INK  •  1,5 
KGN2  •  KC0N2I I  NO  2  I 
Ir  (K0N2)  9077,73,73 

73  KQN 1  -  KC0N1IIND1) 

IF  (KONt)  109,76,74 

74  GO  TO  175,76,77,78,79) , INK 
7  5  TERM  «  UD,’ I'D*  TRIG  1 

GO  TO  90 

76  TEr'M  •  T0E  (  I  TT  l*TR  JG1 
GO  TO  9«S 

7 l  TFRM  «  TCCI ITT|*TRIG2 
GO  TO  90 

78  TFRM  •  TGCI ITTI*TRIG1 
GO  TO  99 

79  TERM  — TH7(ITT)*TRIG2 
99  COUNT  »  COUNT  +  i, 

Sim  «  OlUNnil  #•  TERM 
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mm**"* 


D1IIN01)  -  SUM 
KC3N2IIN02I  -  IB IGI SUM,K0N2I 
IF  (ABSFI TERM/SUMI  -  TOU  101,101,100 
tuO  ic  \ At bt  t rfRM» /TWOKIKON2)  -  •5E-071  101,101,105 

101  K0N1  *  KONI  ♦  1 

IF  <  HSTAN  -  KONI )  102,102,103 

102  KUNl  ■  -N 
GO  TO  108 

103  IF  INSTOP  -  N)  104,104,103 

104  TALLY  ■  TALLY  ♦  l. 

GO  TO  10S 

105  K0N1  >  0 

IF  INSTOP  -  Nl  106,103,103 

106  TAP  »  TAP  ♦  1. 

9077  KONI  »  -1 

103  KC0N1I IN01 I  -  KONI 

109  IND1  •  1ND1  ♦  300 

110  IND2  *  I N02  ♦  900 
INCH  *  I  NO  1  -  1490 
IND2  »■  I  NO 2  -  4470 

IF  IKCOI  129,129,9078 

9078  IFINT1-ITT1  69,69,71 
69  K5T0P»NT1 

7 1  WASTE! KWASTE  I  «  XBARIITT) 

KWASTF-KWASTEU 

124  IFIXWASTE-7)  9124,125,125 

9124  IFIKSTOP-JTT)  129,125,129 

125  XWASTE-1 

IFIKCO)  129,129,126 

126  WRITF  OUTPUT  TAPE  6 , 127 , I  WASTE ( 1 1 ) , I i» 1 ,6 ) 

127  FORMAT  (SX, 6317.7) 

WRITE  OUTPUT  TAPE  6,127,1111  I TH , 12) , I 2« KST ART , KSTOP 1 

WRITE  OUTPUT  TAPE  6,127,(61.1  I  TH,  l  31 . 1 3-KSTART  ,  KSTOP  I 

WRITE  OUTPUT  TAPE  6,177,1C1(  I TH , 14 1 , 1 4- KST ART .KSTOP I 

WRITE  OUTPUT  TAPE  6,127,(011  ITH, 13) ,! 5-KSTART, KSTOP ) 

WRITE  OUTPUT  TAPE  6,127, (HI  I  I TH. 16 1 , 16-KSTART .KSTOP » 

K$TART*K  STARTS 
KSY0P»K$TART*5 
WRITE  OUTPUT  TAPE  6,128 

128  FORMAT  (///» 

129  CONTINUF 

00  135  ICK-l.NTHl 

INCH  •  ICK 
00  135  JCK-1.NT1 

00  134  INK  ■  1,3 
IF  (KCONKINOin  134,140,140 
134  I  NO  1  >  IN01  ♦  300 
133  IN01  •  I  NO!  -  1490 
GO  TO  160 
140  CONTINUF 
C  PARTIAL  SUMS  ON  P 

160  DO  210  !TH«l,NTHl 

00  210  ITT-l.NTl 

INC  «  ITT 
DO  209  INK  -  1, 3 
KON  •  KC0N2(!R0E.TTH,INC) 

IF  (KON)  209,200,200 
200  TERM  *>  OKI  TH,  INC  I 

SUM  «  02(IR0E,ITH,INC)  ♦  TERM 
KON  •  I81GI SUM, KONI 

IF  (A8SF  (TERM/ SUM)  -  TOU  202,202,201 
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201  IF  <  AB  SF I  TERM) /TWOK ( KON I  -  I.F-07)  202,202,2041 
j.ui  vnu  m  «nw  -  u  nuunui 

GP  TO  204 

202  KON  •  KON  ♦  l 

ip  t MCTiP  -  LLOWIKON))  203s 20?: 204 

203  KON  o  -XMAKOF(iLK,iaiGITii6MKbN“JT72>FSUM.On  -  1 

GO  TO  206 

2 04  IF  IKAPP  -  IP)  205,205,206 

205  KON  -  XMAXOF  ( ILK  ,  I BI G!  TERM/S U‘1,0)  ,  2  BIG! TWOK!  K0N~3O72 I/SUM, 01 1 

206  02 ( IROE , [TH, INCI  •  SUN 
KCON2I IPOE ,  I  TH, I NC  )  «  KON 

209  INC  -  INC  ♦  30 

210  CONTINUE 

00  21?  IRCK"! ,NR01 

00  217  ITHCK«l,NTHl 

DO  217  ITCK»1 ,NTl 

INC  -  I TCK 

00  217  INK  *  1,5 

IF  t  KC0N2  <  IRCK, I THCK , INC) 1216,219*219 

216  INC  ■  INC  ♦  30 

217  CONTINUE 
GO  TO  220 

219  CONTINUE 

220  CONTINUE 

221  WRITE  OUTPUT  TAPE  6 ,222 *COUNT,T ALLY* TAP 

222  FORMAT  I  UNITERM  COUNT, F20.C///39H  NO,  OF  N-SUMS  WITH  TOO  FEW  SNAl 
11  TERMS, FI 3.2///34H  NO*  OF  N-SUMS  WITH  NO  SMALL  TERMS, F15.2) 

1220  INC  -  LEAOX  -  l 
LFAP  »  30* INC 

240  00  241  I  ■  1*NTH 

241  THETA! I )  •  57.  2957796*THETA( || 

245  IF  (KOOEN  -  KEY)  247,251,251 
247  INO  -  0 

JUMP  •  >1 
00  249  INK  «  1,5 
00  248  IX  «  l, NT 
INTER  -  INO  ♦  IX 

246  WRITE  TAPE  3, < 102 1 IRO, ITH, I NTER1 ,KC0N2 < IRO, ITH, INTER  I,  IRO  -  1,NR0), 
1ITH  •  l, NTH) 

249  INO  ■  INO  ♦  30 
RFWINO  3 

251  IF  (KODEN  ♦  KEY  -  2)  254,252,3001 

252  00  253  I  -  l, NTH 

253  THETA! I)  •  PHI  ♦  THETA! I t 
JUMP  -  l 

254  CONTINUE 

C  FINAL  OUTPUT 

6001  CONTINUE 

6010  N1  -  NRO 
N2  •  NT 
N3  •  NTH 
LOCI  -  0 

LUC 2  «  12  ♦  LEAOX 

LOC  3  «  3 

INC  1  «  1 

INC 2  •  30 

INC 3  •  3 

LABEL  1  •  0 

LABFL2  •  KFY  -  l  ♦  KOOcN 

6011  INCI  -  INCI  -  N2*INC2 
LOM  «  l 
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661  / 
6  513 

651* 

6315 

6516 

6017 

7019 

6016 

6023 

602* 

6025 

6125 


6926 


6126 

6012 

6013 

6131 
601* 
'  6015 


60*1 


226 

227 

1227 

226 

6115 

60*3 

6016 
6  117 
6118 
6119 
C 


If  I  W  -  II  {  51  2*651**6516 
WP  I  T F  OUTPUT  TAPE  6,6517 
FOP  MAT  I  l  8H 1 0 1  L  A  T  4  T  J  ON  At  WAVF.I 
GO  TO  7019 

WRITE  OUTPUT  TAPE  6,6515 
FORMAT  <  1 IH 1  SHEAR  WAVE  1 
GO  TO  7019 

WRITE  OUTPUT  TAPE  6,6017 

FORMAT  I 46H l SURE  RPOS I T ION  OF  0 1 L  AT  AT IONAL  AND  SHEAP  WAVE*! 
WPITF  OUTPUT  TAPE  6,6019 

FORMAT  <2H  t / t / / / / / // H 6H  THIS  OUTPUT  HAS  THE  FORM  ll\ 

WRITE  OUTPUT  TAPE  6,602* 

FORMAT  (ItH  RHO, XBAR, THETA  III  ///ID 

WPITF  OUTPUT  TAPE  6,6125 

FORMAT  I*3H  THF  SOLUTION  MATRIX  HAS  S RR , S T T  ,  S RT , UR ,UT tt II / 1 1 

00  6119  I l  -  l, Ml 

WRITE  OUTPUT  TAPE  6,6926 

FORMAT  C.MH 

J1  *  i  l  6  LUC1 

WRIT?  OUTPUT  TAPE  6 , 6 1 26  ,  T I TLF ( J 1  ) 

FORMAT  (L15.7J 

IF  <  LABEL  17  6013,6013,6012 

YBAP  «  I  TITLE  I Jl)  -  XSTARI  /  SVMDUM 

WPITF  OUTPUT  TAPE  6,6126, YBAR 

DO  6118  12  "  1,N2 

J2  -  12  ♦  L0C2 

WPITF  OUTPUT  TAPE  6 , 61  31  , T I T LE < J2  1 
FORMAT  (////  f  15 • 7  I 
IF  (LABEL2I  6015,6015,601* 

YBAP  -  (TITLFIJ2I  -  XSTARI  /  SVMDUM 
WRITE  OUTPUT  TAPE  6, 6126, YBAR 
DO  6117  KA  »  1,  N) , 5 
KO  ■  XMIN0FIN3,KA**| 

LIMA  -  KA  ♦  L0C3 
L I  MO  «  KO  ♦  L0C3 

WRITE  OUTPUT  TAPE  6 ,60*1 , ITI TL6 ( J3 ) , J3  •  LIMA,LIMO| 

FORMAT  </////  6X  ,5F21. 7  /// » 

LIMA  «  LOM  ♦  I  NC  3  •  I KA- 1 ) 

LINO  -  LOM  ♦  l NC  3*1 KO-1 ) 

DO  6016  INK  «  l ,5 
IMP  •  KA 

DO  M l 5  LIMP  «  LIMA,LIM0,INC3 

OK  IMP  I  -  02  (L  IMP  I 

KON  ■  KC0N2ILIMP1 

ERROR  *  AB  SF ( T WOK  I KON)  I 

IF  (ABSFIERPOR  -  ,31  -  .5)  226,227,227 

KON  *  -XSIGNFI XINTF(-LOGIOFIERROR) I ,KONI 

GO  TO  2 28  ,,  . 

KON  «  -0 

IF  I  ERROR  1  22e, 1227, 228 
KON  *  NACRCY 
KCONlt IMP!  •  KON 
IMP  ■  IMP  ♦  1 

WRITE  OUTPUT  TAPE  6 ,60*3 , 1  01 ( J 1 , KCONl I J 1 , J  »  KA»K01 
FORMAT  <6X,E17o7,!*,E17»7,I*,El7«7,l*,El7#7,l*,El7#7,l*l 
LIMA  .  LIMA  ♦  900 
LI  MO  ■  LI  MO  ♦  900 
CONTINUE 

LOM  *  LOM  *■  INC2 
LOM  »  LOM  ♦  INCl 
SFLCCTING  MAXIMAL  QUANTITIES 
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280  KWO  «  KWA 

290  IF  (MAX0U1I  300,350,100 

300  WRITE  OUTPUT  TAPE  6,101 

301  (  l  MiiMAAinA  u*tf  A0*K/'/2'in  ab*k  »*LUi»  akc  »ci.u»  mmniri krrrri 

u/n 

OD  310  I  TH*l ,NTHl 

THFT  -  TMETAHT-n 

WRITE  OUTPUT  TAPE  6,302, THFT 

302  FORMAT  ( ////7H  THETA-E14«7  //1H  , 16* ,3HR0E , 17X , 3HSRR, 17X , 3HSTT, 
117X,3H$8  T) 

00  310  IPO-1, NR01 

ROE  •  RHO< IROI 
FMAXt«A2(IPU,ITH,l» 

EHAX2-62I IRQ,JTH,l) 

EMAX3-C2 (IPO, I TH ,11 

I B 1  *  LEAOX 

182  -  LEAOX 

163  -  LEAOX 

00  308  JTT«l,NT 

ITT1  -  ITT 

AMA1-ABSF (E MAX  II 

AHA2-A6SF|A2(IR0,1TH,ITT1 1  I 

IF  <  AHA  2- AKA  1 1  304,304,303 

303  FMAXl«>A2( IRO,ITH,ITTl) 

181  ■  ITT  ♦  LEAOX  -l 

304  AMB1-ABSFIEMAX2I 
AMB2-ABSF<82<  IR0,1TM,ITUM 

IF ( AMB2-AMB l I  306,306,305 

305  EMAX2-B2I IR0*!TH*ITT11 
162  ■  ITT  ♦  LEAOX  -1 

306  AMC l-ABSF IEHAX3I 

AMC2-ABSFIC2 ( IRO,!TM,ITTlll 
IF  <  AMC  2-AMC1  >  308,308,307 

307  EMAX1-C2 I IRO*ITH,lTTll 
IB 3  «  ITT  ♦  LEAOX  -1 

308  CONTINUE 

ETM 1  -  XBARUB1I 
ETM2  -  XBARIIB2I 
ETM3  •  XBARKSi) 

•■•"FT-  rt"TPUT  TAPE  6,300-  ROE  ,  EM  AX  1 ,  EM  AX  2  ,  EMAX  3«FTM1,ETM2,ETM3 

309  I  46 XU*  f  /  i 

310  CONTINUE 

WRITE  OUTPUT  TAPE  6,311 

311  FORMAT  1 16H1NAX IMA  OVER  RHO//28M  RHO  VALUES  ARE  BELOW  MAXIMA/////// 
1/1 

00  319  I TT» 1 ,NTl 

IT  ■  ITT  ♦  LEAOX  -  1 

WRITE  OUTPUT  TAPE  6 ,312 , XBAR 1 1 T » 

312  FORMAT  I////6H  X BAR-E 14. 7// 19X ,5MTHE TA , l 7X , 3HSRR , 17X , 3HSTT , 17X , 3HS 
IRTI 

00  319  ITH«l,NTHl 

THET  -  THE  TAUT*  I 
EMA  X1-A2  (  1 , 1  TH  ,1  TT1 
IB  1  >  1 

EMAX2-B2 1 1 , 1 TH ,1  TT» 

162  ■  1 

FMAX3-C2 (1,ITH,ITTI 
IB3  ■  1 

DO  317  IRO-l.NRO 

IP01  -  IRO 
AMA1-A6SFIEMAX1) 
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AMA2-ABSFIA2  (  IRO. I  ,  ITH,I  TTM 
IFIAM2-AMI)  313,313,3129 
■»l>0  cua  wi  «*2  !  »  on  •  .  t  ru  .  i  rr  j 

IB  l  ■  IPO 

313  AMB l-ABSF ( EMAX2) 

am  a 2  "A B  S  F I  2  2  I !  RO  1 ,  J  TH ,  2  T  T  *  * 

IF ( AMB2- AMB 1 1  315,319,314 

314  EMAX2«B2CIR0l,!TH,ITT) 

IB2  -  IPO 

315  AMC1»ABSF(EMAa3) 

AMC2-ABSFIC2UROI,  ITH, ITT) j 
IF!AMC2'AMCll  317,317,316 

316  EMAX3-C2IIR01, ITH,ITTI 
IB  3  »  IPO 

317  CONTINUE 

FTM1  •  RHOIIBII 
ETM2  »  P  HO  I  I B2 I 
f TM3  -  P HOI  183) 

WRITE  OUTPUT  f APE  6 ,318 , THET , FMAX1 ,EMAX2 ,EMAX3 , ETMl , ETM2, FTM3 

318  FORMAT  (4E20#7/20X,3E20«  7/1 
31°  CONTINUE 

WRITE  OUTPUT  TAPE  6,320 

320  FORMAT  I  1 SHI  MAX  I MA  OVER  THET  A//30H  THETA  VALUES  ARE  BELOW  MAXIMA// 
1//////I 

DO  329  ITT-1.NT1 

IT  -  ITT  ♦  LEAOX  -  l 

WRITE  OUTPUT  TAPE  6 , 32 l , XBAR (I T I 

321  FORMAT  I////6H  XB.  R-E14. 7//17X,3HBH0, 17X, 3MJRR , I 7X, 3HSTT , 17X, 3HSRT 
1) 

00  329  IRQ-1, NR01 

ROE  -  RHOI IPO) 

EMAX1-A2! IPO, 1,1 TT) 

Ifll  -  1 

EMAX2»B2l!R0,l,I TT) 

182  -  l 

EMAX3-C2 1 1 RU, 1 , 1 TT) 

IB  3  ■  1 

00  327  ITH- l, NTH 

ITHl  •  ITH 
AMAl-ABSFI EMAX1 ) 

AMA2-ABSFIA2I IRO , I TH1 , ITT ) ) 

I F ( A  MA  2- A  MA 1 )  323,323,322 

322  EMAX1-A2(IR0.ITM1,ITT) 

IB  1  •  ITH 

323  AMS1*A6SF/IEMAX2) 

AMB2-ABSF IB2I IPO , I TH1 , ITT )  ) 

IF(AMB2'AMB1)  325,325,324 

324  EMAX2-B2 I IRO, I TH l, ITT) 

162  ■  ITH 

325  AMC 1-ABSF1EMAX3) 

AMC2-ABSFIC2I IRO, ITHl.ITT)  I 
IF  t AMC  2- AMC 1 )  327,327,326 

326  EMAX3-C2 I IRO, ITH l, ITT) 

IB  3  •  ITH 

327  CONTfJ'JF 

ETM1  -  THETA! I  81 ) 

FTM2  -  THETA  1 1 B2 ► 

ETM3  -  THETA ( I B3 ) 

WRITE  OUTPUT  TAPE  6 ,32B,ROE , EMAX1 , EMAX2, EMAX 3, ETMl , FTM2, ETM 3 

328  FORMAT  I AE20. 7/2 OX , 3E20. T/ ) 

329  CONTINUE 
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EMAXI-A2I  l  ,  I  ,  1  ! 
161  *  l 
I B  2  -  1 

ii  >  ■  Lr  »0  A 

FMAX2*82<  1,1*11 
1C  1  «  I 
112  •  i 
1C  3  «  LF  AOX 
fMA a3»C2 (1*1*1) 
101  *  1 
102  *  1 


103 

«  LFADX 

DO 

335 

1R01-1 , NROt 

00 

335 

ITHl-1 ,NTHl 

DO 

335 

1 TT1-1 ,NTl 

/.HA  IsABCFlFMAXll 
AMA2>ABSF|A2<iR01,irHl,tTTl)  I 
IF ( AHA2-  AHAl I  Jit, 331,330 

330  FMAXl*A2  f l  RO  l » ITHl* ITT1 I 

mi  ®  if oi 

IK2  -  5TH1 

183  «  ITT1  ♦  LFAO K  -  1 

331  AMB l«ABSP(FMAX21 

AH82-ABSFI62I  I  HOI,  11HI,  ITT  Ml 
IF ( AMU 2- AMB 1  1  333,333,332 

332  FMAX2-B211ROI,  I  TH 1 , 1  TT1 ) 

IC  1  «  IP 01 

IC2  «  ITHl  ' 

1C  3  -  ITT1  ♦  LFADX  -  l 

333  AMC1-AB$F(EHAX3) 

AMC2-ABSFIC2I IP01 , ITHl, ITT tf I 
!F|AHC2-AHCt)  335,335,334 

334  EHAX3-C2IIR01,ITH1,ITTII 

101  •  IR01 

102  ■  ITHl 

103  *  ITTl  ♦  LFAO A  -  1 

335  CONTINUE 

ROFA  ■  R  HO  I  IB  1  I 
HOFfl  -  PHCI  ICll 
ROEC  «  RHOI 101) 

THA  -  THFTA1  ia<l 
THB  ■  THETA (IC2I 
THC  •  THfcTA ( 102 1 
TTA  >  XBARI1B3I 
TTB  •  XBARI1C3) 

TTC  •  XBARl 1031 

WRITE  OUTPUT  TAPE  f ,336,€HAX1 ,ROE A,THA,TTA,EHAX2,R0FB, 

1TH0, TTB, EM AX 3, ROE C, THC, TTC 

336  FORMAT  ( 31H1MAXI MA  OVER  ALL  RHO,THET A,XBAR// ft 

15H  SPR ■£ 14, 7,6X, SHF OR  RH0-F14, 7,6X ,6HTHETA»E 14,7,6X, 5HXBAR «F14 , 7// 
25H  STT«F14.7»t  X, SHF OR  RH0-E14. 7 ,6X .6HTHET A«E 14. 7 ,6X, 5HKBAR-E  l<>. 7ft 
36H  SPT»E 14*  7,6X «  RHFOR  RHO- t 14. 7 , 6X , 6HTHETA»E 14, 7 ,6 X , 5HXBAR «F1 4, 7/ I 

350  IF  ( JUMP  I  508,1001,351 

351  JUMP  «  0 

INC  »  NX  -  NT 
00  364  ICX  ■  1,NT 
IX  «  NT  ♦  1  -  ICX 
00  3 5A  IRO  «  l , NRO 
INf)  •  30*1  X  -  30  ♦  IRO 
00  353  ITH  ■  l.NTH 
00  352  INK  *  1,5 
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IN'ER  »  INO  *  LEAP 
C2IINTLRI  ■  02  (I  NOI*AMPFAC 
KCON24 INTFR |  -  K COM2 (INO) 

352  INO  ■  INO  ♦  90 0 

in  i  tu n  -  lun  -  /./.a- 

354  CONTINUE 

IF  UNCI  3001*341,357 
357  DO  360  IX  -  1 , t NC 
00  360  IRO  •  1  * NRO 
INO  »  30* 1 X  -  30  ♦  IRO 
00  350  ITH  *  1* NTH 
00  354  INK  *  1,5 
D2UN0I  *  0. 

KC0N2< INO I  -  -1 
356  INO  «  INO  ♦  900 
354  INO  «  INO  -  4447 

360  CONTINUE 

361  FIX  •  -l, 

00  370  INK  «  1,5 

If  ( INK  -  21  362,363,362 

362  FIX  *  -FIX 

363  DO  370  I X  *  l  ,NX 

RFAO  TAPE  3  « (I f I IRO, ITH) , JF ( IRO, 1THI  , IRO- 1 , NRO) , ITM«1 .NTHOLD » 
IF  INF  NTH  I  3001,565,364 

364  DO  365  ITh  •  1,N*NTH 
JTH  -KE VTH I  I TH ) 

INDEX  -  ITH  ♦  NTHOLD 
00  364  IRO  ■  l, NRO 

F  I  |R0,  INDEX)  -  ’r  I  X*F  1 1 RO ,  JTHI 

365  JFI IRO, INDEX)  =  Jfl IRO, JTHI 
565  INOEX  -  IX  ♦  30*1  INK-1 1 

DO  360  ITM  •  1 , NTH 

DO  366  IRQ  «  l, NRO 

00  -  024 IRO, ITH, ( NOE XI 

DF  -  FI  IRO,  ITH  I 

KON  «  KCON24 IRO, ITH, INDEX) 

JON  ■  JF I IRO, ITH I 
LON  -  -l 

If  I IkON-XABSF (KON)  I  *  I JON-XABSF ( JON)  1 1  367,366,367 

366  LON  ■  1 

367  KON  -  XA0SF4KONI 
JON  -  XA0SF4JON) 

A0  *  MAX1F(A9$Fi  DC* T WOK  I KON I  I , A9SF(0F*TWQK( JON) I  3 
00  »  DO  ♦  DF 

KC0N2IIR0, ITH, INDEX!  «  LON  *  ISIC4 A0/DD,OI 
360  021 IRO, ITH, INDEX)  -  DO 
370  CONTINUE 
NT  «  NX 
LEADX  •  1 
KODEN  -  2 
DO  TO  6001 
C  ERROR  RETURNS 

3001  WPITE  OUTPUT  TAPE  6,3002,ENS,KAPN,IP,!RUE 

3002  FORMA 1  4 15H1  ERROR  RETURN  ft  «H  ENG»E14#7  /  6H  K4PN«I6  / 

14H  |P-I6  /  6H  (ROE-16  I 

3003  CALL  POUMP 
GO  TO  1001 

3006  WRITE  OUTPUT  TAPE  6,3007,LEROR 

3007  FORMAT  419HIERROR  IN  INVERSION  It  7H  LER0R-I6) 

GO  TO  3001 


SUBROUTINE  ZflHSP  ( P, E2 , 65 , E6 , FI  1 , 0 3T IL ,Q3BAR ,EPS6, FPST , VMOUM 

I  .K [T  VPC  jtr  APP  .  «nn<  0  .r  »  Vf  N  .  Tl  THT  .  TlTT  .  T  I  T  .  T9  f  5  OCJ  .  CM  I  .  CU.  *r  Apr  » 

DIMENSION  Pt  J( t)^EHj(iliEM<12fli 
FXPON-EXPF I-03TI L ) 

TP  IGl«<INF<t)3HAR  t 
TP IG2-COSF { Q3BAR  ) 

4R-fPST4lftEJmTRF.il  IHHI.OtFPSB)  *Et4J(  l l-li,0-EPSH)*EMJI3) 
AK  *0.  4  *A  R 

AI«EPST*IEMJl IItFMJI 3)1-1 1.0tEPSB)*REJI 11*1 l.O-EPSBI*REJ5 31 
Ai-0.**AI 

T!LDurf«F  6/ IF  2*  P*  3,141592  7) 

ATfi  ■(*  TIL0UM>*<Q3BAR*EMJI2)-03TlL*REJI2>t 
ATI-Til.nUN*l03BAR*  £J(2)tQ3TII*EMJI2)  ) 
f  Ml  1,1  1-0.0 

FM(2,1 )«FXP0N*<A«*TRIG2*AI*TR!Gl ) 

EMI  3, l)« l-FXPONI  R|ATR*TRIG2+AT| *TRIGl) 
f*H4,l).0.0 

EM ( 4 , l ) «  I  EXPON) 4(AI*TP|G2-AR*TRIGl 1 
FM| 6, l )  -l-EXPONI  <MAT!*TRIG2-ATR*TR|G1) 

PI-3.1415927 

pp |>P*PI 

IMKODSP)  1001.1001,1002 

1001  SP-i.O 

GO  TO  1003 

1002  CAPP-KAPP  ♦  10 
PPIP-PPJ /CAPP 
SP«SINF<PPIP)/PPIP 

1003  PPI2-PPI/2.0 

GO  TO  (900l,9002r9003l,<TyPf 

9001  AP-I  2..  O/PP  I  ?  )  *51  NF  |PPI2  )  *S  INF  I  PPI2*Fl  1 )  *  SP 
AP-42,0/PP12)*S!NF(PPI2)*SINFlPPI2*Ell)  •  S® 

GO  TO  9004 

9002  AP-I (C AVlN4l«0)FCAYINI*$INF«PPI*T|T)»5INF{PPI*T0T) 
l-l 1.0/CAYIN)*SINF(PPI*T2T) 

AP-(2.0/IPP1**2*T1TOTMAAP 
AP  ■  AP  *  SP 
GO  TO  9004 

9003  C0M1-SINF(PPI*T1 Tl 
C0M2-C0$F(PP|*TIT) 
fcXPl-FXPFlACAPT*!  1.-T1T) > 

API-I2./IPPI**2*T1T0TI)MC0HI“SINF(PPI*T0T ) ) 

AP2-2, /PP1 
AP3-ACAPT442tPPI 4*2 
AP4-AP2*C0M2 
IP-P+.05 

EXP7-1 (• l.01**!P»4PPI 
AP5-I2./I AP34EXP1 II *EXP2 
APfc*l2»/AP3)4i AC  APT*C0MItPPI 4C0M2I 
AP-API-AP4-AP5TAP4 
AP  «,  A/  *  SP 

9004  UF-I 1.0/PPI I*<E6/IE5*VMD0M4E2) ) 

00  9005  1-1,6 

9005  EMI  1 , 1  )-F M|  I  ,  1  )*  AF 
EM|3,1J.UF*EMI3,1. 

FMIE.t )-UF*6M<6, ll 
RETURN 

ENO 
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SUBROUTINE  ZKHS*  (P,E2.E5.E6,Ell.04TH..04BAR.FPSB.FP<;T.»/*ftn« 
l.KTVPf ,MPP,KOOSP.CAyiN,TlTOT,TOT,TlT,T2T,RE J, EMJ, FM, ACAPT I 
DIMENSION  RE J< I) 

EXPON  ■  EXPEI-04TIU 
TRIG1  *  SINEIQ4DAR) 

TRIC2  *  COSE  I04BAR ( 

CR  *  -EHJI3I 
Cl  •  REJOl 

AOJ  -«  .225  079074  *  £6  /IP  ♦  E2  *  SQRTF I VMDUM  1 1 
CTR  »  AO  J  *  I04T I L  *RE J I  2 t  -  O40AR*ENJ( 2* » 

C  T  I  *  AO  J  *  IQ40  AR*RE  J 1 2  I  ♦  Q4TIL*EHJ<2M 

E«m  «  o. 

EH ( 2 1  *  EXPON  *  (CRPTRIG2  CI*TRlGi) 

EM<3>  *  -EXPON  *  (CTK*TRIG2  ♦  CT I *TR 1GI I 
f  H  (  4 1  «  0 » 

EM(4 1  *»£ XPON  *  (CRPTRIGl  -  CI*TR!G2> 

EM  ( 6 )  -  FXPCIN  *  (CTRPTRIGI  -  CTI*TRIG2» 

PPI  •  3.  14  1  5  92E5*p 

IE ( KOD  SP 1  100' ,1001, 1002 

1001  SP-l'O 

GO  ro  1003 

1002  C A  PP  «K APP  ♦  10 
PPIPoPPI/CAPP 
SP«SINflPP»P)/PPIP 

1004  RPf  ?«»P  J  /?.  0 

GO  TO  19001,901)2.90031  ,KTVPE 

9001  AP*I2«0/PPI2)*SINF(PPI2)*SINE1PPI2*E111  *  SP 
GU  TO  9004 

9002  AP-I  (CATIN*l.OI/CAVINI*SINF<PPI*Tin-SINFIPPlPTOTI 
1~(1.0/CAV1N)*SINF1PPI*T2T} 

AP«I2,0/<PPI**2*T1T0TI )*AP 
AP  •  AP  *  SP 
GO  TO  9004 

900)  C0H1«SINE1PP1*UT) 

C0H2»C0SFIPPIPT1TI 
FXPt»FXPF IACAPTPI l.-TITI) 

AP1-I2./(PPI<»P2PT1T0TI  )*(C0M1-:!N(  !~?I*T(r  ♦) 

AP2-2./EPI  -  .  _ 

AP3»ACAPT**2aPPI **2 
AP4»AP2*C0H2 
IP*P*. 05 

EXP2-I I-1.0I**IPJ*PP! 

AP5*12./IAP3*f XP1»I*EXP2 
A''4.I2./£P3I*(ACA.,»T*CONUPPI*COH2I 
AP-AP1-AP4-AP5+APE 
AP  •  AP  *  SP 

9004  l)E  ■  16  /I  PP  12  *  E2  *  E5  *  $QR TE( 2. •VMUUM)) 

on  9005  1*1.6 

9005  E M (  1  I  ■  FMm  *  AP 
EM  1 31  »  EMI3)  *  UF 

EHltl  m  EHI6)  *  UF 

RETURN 

END 
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SUBROUTINE  RHSP  IP^.FS.EE.Fli.l^TIUQaBAR.FPSBiFPST.VHDUM 

I  >VT  VOf-  .MAOO.vnftCO.r  4  v  f  fu  .  ▼  1  fdT  T  nr  yit  r^r  .»  r  *  r  w  .  «n  a*  «.  ^ 

-.  ---  .  •  *  -  —  —  -  *  ---•-»  •  -  ■  -  •  »  ■  »  ■  r  ■  •  ■  »•»-<  *  ■'  >.  w  t  •  ri  V  ?  »»vm  r  i  f  |tf  r 

0  IMFNS  ION  PF  J<  1)  ,E4jm  ,EM<12«li 
NFXP»( N»l»/2 

B*£  TOD  <•«  -  I  .  /)}  «r*MC  wp 

PlnV2-li 570796327 
P  I  *3,  1 41  5927 
1F(N-2*(N/2II  1,1,2 

1  Bl-REJiN+l) 

82*FM J  < N-f  1 1 
Cl— FMJ1N+21 
C2  — FMJ( N( 

C3»PCJ(N*2) 

C4«PFJ  (M 

01«RFJ(N*31 

02«BFJ1N-1» 

03«F  M J { N* 3  I 
OA«rMj(N-l> 

GO  TO  5 

2  di*E«j(im> 

02— «E  JIN*|» 

C1*RFJ1N»2» 

C2-RFJ1N1 

(.  3«ENJ  <N*2  I 
C*-E0J(N) 

01-  FMJ1N*3» 

03  — PF  J1N*31 
1FCN-1I  3,3,4 

3  02  —  FHJI21 

D4«*pe JI2I 

GO  TO  5 

4  D2w*<=MJlN-ll 
04— RF  JIN- 11 

5  f l«f  AC  TOM*P  I 
F2»FACT0R*P10V2 
A8N-FI*B1 

ATM»F1#«2  -  - 

ABCNi»E2*lCi-C2) 

ATCNI-F2*1C3-C4> 

AHSNI-1-F2»*1C1*C2I 

ATSNI"1-F2)*1C3*C4| 

AeCN2«l~F2l*COI*D2l 
ATCN2«1-F2»*(D3*04I 
A6SM2»F2* ID  1-021 
ATSN2»F2*103-04i 
F3»0.31«3j9B9 

APNR.gPSf*! ABN-49CN21 »l l.O*FPSB)*ATNM l.0-FPS9)*ATCN2 
APNR  »F 3* APNR 

APNI«FPST*l*TN-*TCN2»'U.(J*EFS0»*AdN-l  1.0-£PS9)*ABCN2 
APN I *F 3*APNI 

CPNP-F  3*4FPST*ABSN2-(l,0-EPSBI*ATSN21 
CPNI*F  J*IEPSt*AISW2*l 1. 0-EPS9)*A9SN2» 

P4.|2.0*Ff )/lP*PI*PI*E2» 
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ATPNR*F6*(  G3BAR*ABCNl  *03T I L*ATCN1  I 
ATPN J  «F**{ Q3BAP*ATCNl-33TI L*ABCNl i 
CTPNR«( -F6>*(O30AR*AaSNl*Q3T!LPATSNl » 

C  f  rni-  i-r  •» »  »»  jnil-W  »i  1L»  ABir,i  I 

FXP0N*EXPF(-03TI  LI 

TP JG1-SINFI0SBAR) 

in  iG2*ClJ6M  WABAK  I 

FM<1, 11*0,0 

EM(2, 11*0.0 

FHI  3,ll*f XPON*  (APNR*TRIG2*APNl*TRIGll 
EM(*,ll«EXP0N*(CPNr*TRIG2*CPNI*TRIGl» 

EM ( 5 , 1 )  ■ ( -f  XPONI  *<  ATPNR*TRlG2«-AfPNI*TRlGl  I 
EM (6,1  l*(-EXPON)  *(CTPNR*TRIG2«-CYPNIPTRIGi; 

EM(7, 11*0.0 

em(b, n-o.o 

EM (9,1 l«l  F XPONI *<APNI *TR1 G2-APNR *TR IG l I 
EM  ( 1C,  11*1  EXPONI*(CPNJ*TR!C,2-CPNR*TRIGll 
EM(U,  ll«(-EXPONI*(ATPNI*TRIG2-ATPNR*TRIGi» 

EM<  12,  1)«(-EXP0NI*(CTPN1*TRIG2-CTPNR*TRIG1» 

PJ-3.1415V27 

PPI*P*PI 

IF(KODSP)  1001,1001,1002 

1001  SP-1.0 

GO  TO  1003 

1002  CAPP-KAPP  ♦  10 
pp|p*ppi /CAPP 
SP«$INHPPIP)/PP1P 

1003  PPI2«PP1  (2* 0 

GO  TO  (9001, 9002, 90Q31,XTYPE 

9001  AP *<  2, 0/PPI2 l*$I NF( P PJ2 1 *$INF( PPJ2*E 11 1  •  SP 
GO  TO  9006 

9002  AP-({CAYIN»l.O)/CAYiN>*(|NF(PPI*TiTI-$INF(PPI*TOT) 
1-<1.0/CAYINI*SINF  *MI<M2TI 

AP-(2.  0/(PPl*<»2*rtT0TJ»*AP 
AP  r  AP  •  SP 
CO  TO  900* 

900  >  .rJMl«SINF(PPl*TlTI 
CC’N2*C0SF(FPI*Tl  T  > 

EXP1-E XPF( ACAPTP( 1.-T1TI t 

API*(2,/(PPI**2*TITQTJJ*| CQM1-SINF IPPJ  *T OT ) > 

AP2-2. /PPI 

APS* AC APT**2*PPl **2 

AP**AP2»C0M2 

|P*P*.  05 

1  XP2*( (“1,01 •* IP! *PPI 

AP5*<2./(APS*E XP  1 1 1 *EXP2 

APf  *(2./AP3J*(*CAPT*C0MUPPf  *C0M2  I 

AP-AP1- AP6-AP5  *APf 

AP-AP1-AP*-AP5*AP6-AP7 

AP  •  AP  *  SP 

900*  OF *| l.O/PPJ l*(€6/(E5*VM0UM*E2t) 

00  9005  (-1,12 

4005  E«<I,lt-EM(!,l(*AP 
€M(S,ll*t>F*EM<5,  II 

€MU,1I-UF*€M<  6,11 

€M'll,l.*Of*EM(ll,l» 
fM112, l)"0f*EM(l2,ll 
RCTURN 
f  NO 
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SUBROUTINE  RHSH  (  P  .  E2  .  E  5  .  E/> .  F  1  1 . 04  T  1 1  .  04  H»o  .  to  <:o  .  i-  d  c  r  .  vmaijm 
1,KTYPF,KAPP,K0DSP,CAyIN,  T1  TOT, TOT, TIT  ,T2T,RPJtEMj,"ACAPT,NrEMI 
DIMENSION  P  E  J I  1)  ,  EM J ( 1  ),EM< 12,1) 

NFXPMN*; 1/2 
PAC.TOP*{~1»OI**NEXP 
PiOV/-»l,  570796327 
PI*3, 1415927 
IE (  N-2*( N/2 ) 1  1,1,2 

1  Cl  — EMJ(N*2) 

C2— EMJ< N I 

C  3*Rfc ) IN+2 ) 

C,4*P€  J  (N ) 

PI  -»RE  JIN+JI 
02  =RE  J ( N~1 1 
D  3*F  M J I N*  3 1 
D4-FMJ ( N— 1  I 
GO  TO  5 

2  C  l  BEJ(N*2) 

C  2*F.E  J  ( N  I 

C  3»f  MJ  ( N  +  2 ) 

C/,bC.MJ  (N  I 
01«  EMJIN+3) 

03— RCJ(NO) 

IEIN-1)  3.3,4 

3  02  — EMJ<2) 

04— REJI2I 
GO  TO  5 

4  02— FMJIN-ll 
O'— RFJIN-1I 

5  F 1 »F  AC  TOR*P I 
F2-FAC TUR*PI 0V2 
A8CNI«F2*(C1~C2I 
ATCNl«F2*(C3-C4l 
A8SNl-<-F2|*<CUC21 
ATSNl«(-F2|*(C3tC4J 
ABCN2“<*-F2l*(Dl«-Q2) 

4TCN2*<~F2|*<D5M>4> 

A6SN2-F2*<01-02I 
AT$N2»F2*< 03-041 

F 3  •  . 6366198 
APNR  «  F3  *ATSN2 
APNI  «  -F3  *A8  SN2 
CPNR  «  F3  *ATCN2 
CP  MI  «  -F3  *  6BCN2 

F4  -  • 143209793  *  E6  /IP  *  E2  *  SORTF< VMDUM) ) 

ATPNP  -  F4  *  (04BAR*A8SN1  ♦  04T I L* ATSN 1 1 
ATPNI  «  F4  *  <04BAR*ATSN1  -  09  TI L*ABSN 1 ) 

CTPNR  *  ( 04 8 ARP ABC N1  »  Q4T I l*ATCNl I 

CTP.NI  «  F4  *  I046AR4ATCN1  -  a4T)^*A8CNlI 
FXPON  •  tKPF(-04TIiU 
TRIG1  «  SINFIU4BARI 
TR IG2  -  C0SFI049 AR) 

FMI 1,1 )»0.  0 
EMI2,1I«0,0 

F*M  3, il-EXPON* <APNR*TR1G2*APN|4TR1G1 I 
EM  <4,1 1 -E  K'PON*  (CPNR*TRIG2+CPNI*TRIGi  I 
F M (  5 , 1  I ■  (  *E XPGN)  *UTPNR*TRIG2+ATPNI*TRIGU 
EMI  ‘..,lI«t-F9  10N)*I.CTPNR*TRIG2+CTPNI*TR  IG1  ) 
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F  M  I  7 . 1  I «  0.  n 
EM<8, 11*0.0 

FM  (  9, 1 1  *  (  EXPON)  *<APNl*TfltG2-APNR*TRIGl> 

FMi  I  o.  1 1  ■(  f  xponi  t  *«rPM i  *tb  ;  r-2-CPNP*TP !  Gi  ! 

EMC 1 1 « [) *C -EXPON I *< ATPN I *TRI G2-ATPNR*T  R I Gl ( 

EMC  12 1 1 I *< -EXPON )*CCTPNl*TRIGC-CTPNR*TR IGl  I 

PPI.P*PI 

tF(KODSP)  1001,1001,1002 

1001  SP-1.0 

GO  TO  1003 

1002  CAPP«KAPP  ♦  10 
PPIP-PPI/CAPP 
SP-SINFCPPIPI/PPIP 

1003  PPI2-PPI  /2.0 

GO  TO  <9001, 9002, 90031, KTYPE 

9001  AP*C2.0/PPI2l*SINF<PPm*S!NF<PPi2*Ell)  •  SP 
GO  TO  9004 

9002  AP»HCAVIN-tUO)/C£YIN)*SIttF(PPI*TlT  l-SINFC  PP  I*TOT  I 
1-C l.0/CAYINI*SINFCPPI*T2T> 

AP”C2c  0/CPPI**2*TlT0TI l*AP  ' 

AP  «  AP  *  SP 
GO  TO  9004 

9003  C0M1*SINFCPPI*T1T| 

C0M2*C0SFC  PP  f  *Tl T) 

EXPl-EXPFCACAPT*  Cl.-TITI  I 

AP1*<  2./ (PPi**2*TlT0T) I  * ( C0M1-S I NF ( PPI *TOT  M 

AP2«2./PPI 

AP3-ACAPT**2*PP[**2 

AP a*AP2*C0M2 

IP-P*.0f 

EXP2«((-l.0)**lP»*PPI 

AP?«<2./CAP3*FXPll)*EXP2 

APfr*( 2./AP3  )*<  AC  APT*C0M1*PPI *C0M2 ) 

AP "API-AP4- A P5tAP6 
AP  ■  AP  *  SP 

9004  UF  ■  E  6  /<  PPI2  *  F2  *  E3  *  SQPTF <2,*VMDUM)| 

DO  9005  1*1,12 

9005  EMC  I,  l l-EMC  I  ,1  l*AP 
EMC5,n-UF*EM(5,  II 
EMC  4,1 l*UF*EM<6, II 
EMC  11,11 *UP*EM ( 11,11 
EM! 12,11 «UF*EM( 12,11 
PF  TURN 

END 
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SUBROUTINE  LINSYS  ( A i Y  *  H, L ! 

01  HF NS  l  ON  A  ( 12  *1  21  ,YU21 
HI  ■  M  -  l 
00  150  K  »  l,Ml 
KP  *  K  ♦  1 
X  -  0, 

00  110  l  ■  K» M 

IF  (X  -  AB  SF I A  1 1 ,KI I  1  100,110, liO 
100  X  *  AB5F<Ai!,K)> 

L  -  I 

110  CONTINUE 

IF  (X)  120,120,110 
120  1*0 

GO  TO  190 

130  00  1A0  J  -  I ,  H 
X  «  A  (  K ,  J I 
A I  K  *  J )  «  AIL, J) 

140  A(L,J)  -  X 
X  «  YIKI 
Y(K)  •  YIU 
YIU  ■  X 

DO  150  1  »  KP  »H 

X  •  A ( 1 , K )  /  AIK, K> 

veil  ■  vm  -  vi ki *x 
00  150  J  -  KP,H 
150  A I  I , J I  -  All  ,  J  )  -  A( K, Jt  *X 
IF  I  A IN, HI )  160,120,160 
160  Y I  H )  «  Y  I  M  I  /  A<  H,  M  ) 

K  •  H 

00  180  I  ■  1,M1 
X  •  0, 

KP  ■  K 
X  »  K  -  1 
00  170  J  ■  KP ,H 
170  X  ■  X  ♦  A I K , J ) * Y I J I 
180  Y1K»  ■  I VI K I  -  XI  /  AIK.KI 
190  RF  TURN 
ENO 
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SUBROUTINE  TILDE  < K, 0,R OE , THETA, RE JT IL , EHJTI L, 

IRE  VT I L ,F  MVTI L»ENG ) 

01  MENS  ION  REJTIU1I  .EMJTIL(l)»REVTlL(l|fEMVTIL(l> 
IMORE-1 


tc  i  rue  n  i  ii  ^ 


1  IMORE-2 

2  N  •  ROE 

tc  Itnn  rturnti  «i  •  <«mmi  «a«i 

•  •  AwvwriwbM  t  ua  UHnrtbl  IWUI  UUtM 

2001  N  •  10  ♦  XMA XOF ( K » Nt 2 0 1 
810  RE  JTIL IN+2  M0*  0 
EMJT1LIN42I-0.0 
REJTILIN+ll  •  l.E-18 
EMJTlUN+ll  -  0. 

2003  RO-ROERQ 
Fl-1.0 

GO  TO  (3,4) , MORE 

3  X«RQRC0SF|THETAI 
Y-RQRSINFITHETAI 
FI-X/P0 
F2«V/|RQ*ROEI 

4  Fl-Fl/ROE 
QSQA»C**2/4.0 
00  if  I » 1 »  N 
L»Nt1-I 

EL-l 

REJTILIL)«F1RREJTILIL«1»-IQS94AREJTILILr2)  >/ (EL*t ELRl.01  I 
FMJTILILI  -  0. 

GO  TO  15, 61, MORE 

5  REJT!lIU»REJTIUl)*F2*ENJTIL(L*l> 

FMJTILI L  >«F 1REMJTI LI L+l l“F2*RE JT IL I L*1 t- I0SQ4*EMJT IL I L+21 1 
1/ ICL*  t  ELT 1*  OH 

6  IF  lABSFIREjTILlLM  ♦  ABSF  I ENJT IL (  L  1  I  •  1.E30I  16,16,2004 

2004  00  2005  J  -  L,N 

REJTILIJ)  ■  REJTIHJI  R  UE-30 

2005  ENJTI L ( J I  -  EMJTIHJ)  •  1. E-30 
16  CONTINUE 

ENG"0*0 

CALL  JTILO  IQ, ROE, THETA, FP,GP,ENGI 
IF1ENGI  7,7,28 

7  GO  TO  (8, 4|,  MORE 

8  A  ■  RfJTIUn 
B  -  FMJTIU1) 

HIGH  -  MAX  IF  I  ABSF  ( Ai  ,ABSFI  BH 

AH  ■  A /HIGH 

BH  •  B/HIGH 

OENOM  •  A»AH  ♦  B*6H 

A  ■  IF PR AH  ♦  GPRBHI  /  OENOM 

B  -  IGPRAH  -  F PR SHI  /  OENOM 

GO  TO  10 

9  A-FP/REJTILI II 

10  00  12  I ■ 1 ,N 

DUHl-REJTILlIt 
RF  JT  IL  1 1  MAROUNl 

GO  TO  1 11,12},  MORE 

11  OUM2*FNjTiL(I> 

REJTIHI  I-REJTILIIM8R0UM2 
FMJTILI |  )*A*0UM2 TB*DUM1 

12  CONTINUE 

CAP  Aa-0, i 1593131 57+LOGF ( RQ) 

SUMRE* 10*63661 97721 R(CAPARFP>THETARGP| 

SUM  M"  10*  63661 *7 721 • (C APARGP+TMET ARFP) 

GAMrI. 273239544 
F0UMMY-QS04 


FACTOR «FOUHMY/2. 0 
MDV7-I N« I  I /3 
KTEST-l 
LTFST-1 

nn  »4«i.wny2 

M2-2*M*1 
TERM»GAM*P ACTOR 
GO  TO  <15,1701 J.KTEST 
15  TERHRC»TFRM«RFJTIL(M2I 

IF(A8SF(TERMRFI“  UOE-ll)  46,46,17 
Aft  KTFST-2 
GO  TO  1701 

17  SUMRF«SUHR£+T£RMRf 

1701  GO  TO  ( 1702  ,22 I » LTE  ST 

1702  GO  TO  ( 19,18) « IMORE 
IS  LTEST-2 

GO  TO  22 

19  TERMlM«TERM*f MJTILIM2) 

If <ABSF(TERM!M)-1.0E-ll»  20,2 Of  21 

20  LTfST-2 
GO  TO  22 

21  SUMIM«$U*IM*TFRMIM 

22  IFUTESTHTE  ST-41  23,25,25 

23  6*-<M2*l)*M2 

FAC  TOR •( FAC  TOR  *F  OUMMV I  /BM 
EM»M+l 

GAH-(-GAM)*((EM-l,0»/EM| 

24  CONTINUE 
€NG“3«  0 
GO  TO  28 

25  REYTIL11 )»SUMRE 
EMYTILI 1  •■SURI N 

GO  TO  (27,261, IMORE 

26  T0R*QSQ4*SUNRE*REJTIL(2  )-Q/( 3a 141 59265 3*RQ) 
REYTlim-TOR/REJTIUl» 

FMYTIL12I  «  0. 

GO  TO  2009 

27  BOTTOM wREJT IL<  l)**2*EMJTILm»*2 
Ti«RejTIL(ll*RFjmC2KEMJTIUU*ENJT!L(2» 
T2*RE4TIC(ll*EHJTILI2  )-EMJT  IL(l)*R£JTIL(2t 
T3-1.0/(3, 14159265*R0E*R0I 

T4«x*r  ejtil  ( 1  )-v*emjti  un 

T5“Y*R  E JTIL ( 1 ) ♦X*EMUTI L ( 1 1 

REVTlL(2»«(0S04*R6VTIUl»*Tl-0Sa44EMVTIL(l|4T2 
1-T 3* T«  > /BOTTOM 

EM YT IL  <  2  l«(t)$Q4*R6  YT IL  ( 1 1 •T24QSQ4*ENYT tl  ( 1 ) *T| 
1+T3RT5) /BOTTOM 

2009  IF  ACCUMULATOR  OVERFLOW  2010,28 

2010  ENG  •  -2. 

28  RETURN 
END 


SUBROUTINE  BESSEL  1 K, Q, THETA, RE J , EMU, EN3I 
DIMENSION  REJIll.fcMJIl) 

10  FORMAT  (24H  OVERFLOW  IN  BESSEL.  0  -EU.7,9M  THETA  *fcl5.7,TH  KAPN 
1-131 

If  ACCUMULATOR  OVERFLOW  100 >100 
*  GO  K7E  571  *■  Q  ♦*  *G» 

KTEST2-K 

JFlxTF.  STI-KTEST21  1,1,2 

1  N-KTE  ST2- 10 
00  TO  3 

2  N*KTESTltlO 

A  REJIN-21 *0,0 
EMJ4N-21-0.0 
REJIN-ll-l.UE-37 
EMJINMI  *0.0 

X»0*C0SF  4THCTA  I  — 

V-0*  SI NF ( THF  TA  1 
Nl-N+1 

osc-c*o 

01-12. D*XI/OSO 
02-12. O-VI/OSQ 
00  4  I « l ,  N 

L-Nl-I 
EL-l 

RLJILl-tL»IOl*RtJIL*l)AQ2,*EMJIL*ll  l-REJ(L-2» 

EMJILI -EL*IOl*  EM  JIL -11-02*86  JILT  l 1 I-EMJ1L-2 ) 

IF  I ABSF  t  RL  J  I L  )  I  ♦  ABSF I EMJi  LI  I  -  1.06301  4,13,13 

11  00  14  J  -  L,N 

RE  MJ)  •  REJIJI  *  1.06-30 
14  EMJIJ1  •  6MJIJ1  *  1.0E-30 
4  CONTINUE 
ENG-0.0 

C  CALCULATE  ALPHA 
KF-1 
KG- 1 

R  TEST-0.  00000003 

TP-REJIJI 

CP-EMJI3I 

DU  *007  I  SUM  -  3,N ,2 
00  TO  19001 ,*0031 ,KF 
4001  PP-FP-BE J( I  SUM) 

RELF-ABSFIAEJI  ISUMI/FPi 
IFIR6LF-RTLSTI  9002,9002,9003 
9002  KF-2 

9001  CO  TO  19004 ,9006 I, KG 

9004  0P-6P-6MJ1 1  SUM  I 
RELC-ABSF1EMJI ISUMl/OPI 

IF  IRELG-RTESTI  9005,9003,9006 

9005  KO-2 

9006  IFIRC4KF-4I  90  0  7,9006,9008 

900  7  CONTINUE 

ENG  -  -2. 

60  TO  200 

9008  A  -REJ1  11*2. 0*FP 
H  -EMJI 1 l»2.0*GP 

3005  HIGH  -  MAX IF  I ABSF  4  At , ABSF I H) I 
X  ■  A/HIGH 
V  -  B/HIGH 
OENOM  -  A*X  ♦  tt*V 
A  ■  X/OENOM 
U  ■  V/DENUM 
4006  00  e>  I  -1  ,N 
DUMI-RE J( 1 1 
0UM2-EMJ  ( 1 1 
Rt  J(  1 1 -A*UUMI*B«  l)UM2 
LMJt 1 1 -A*0UM2-B*0UM1 

6  CONTINUE 

IP  ACCUMULATOR  UVCMFLQW  200,7  . 

200  HPITf  OUTPUT  TAPE  U, 10,0, THETA, K 
ENG  •  -1. 

7  RETURN 


saaSftaXi**  —  - 


•U*.  -■  .k 


su^huui  iNfc  J  1 1  Lu  <QP  fKOtP»THETAP,FP, 
DIHFNtfiiN  a  (  M  .  nu  i  i  t.  Tut  r  a  <  i  ,  ,t ,  , 

-•-»»  M»>  I  A  t  i  f  f|  V  J| 

Of  i)  *UP 

ROFOI-FGEP 

T M fe  T  A  (  1 1  «TH£  TAP 

0 ( 2l«U.O 

RQt<2>«0.0 

|MPTAi2)«0.0 

KTf $T«l 

LTEST'l  - 

IF (THETAJ  2,1.2 

1  LUST-2 

2  F*  1.0 
C*0.  0 
A*- 1.0 
F«Q*R(>£/2.0 

00  5  I  >1,1000 

P-1.0 

FI -1 

00  J  J-I.l 

rj-j 

P-P*(«/tJI*(B/tJI 

3  CONTINUE 

GO  TO  001.3041,  KTEST 

301  PF-P*CUSH  2.0*FI  *TrtETA| 

IF  I ABSF ( PF 1-0. 1**1 11  302, 302,303 

302  KTEST-2 

00  TO  304 

303  F*F***PF 

304  GU  TO  1303,308) ,  LTCST 

305  PG*P*SINF  <2.0*EI*THETA1 

1F<  A9SF  (PG)-O.  l**U»  306,306,307 

306  LTEST-2 
GU  TO  306 

30 T  G-u*A*PG 
308  A*-A 

IF|kTE$T*UEST-4I  5,503,603 
5  CONTINUE 

GO  TO  <501, 5021, KTUST 

501  ERROR-l.O 
GO  TO  504 

502  £PB0R«2.U 
GO  TO  504 

503  FP>F ( l 1 
GP«G< II 

504  RETURN 
END 


.  xV. ' 


-riA  "  . 

•;*  4- .  .j&k  mk4r*m>. 


FUNCTION  ISIGM'Nt 
IF  U»  100,140,100 
100  V  ■  1. 44269504  •  S.OOFI  ABSFIA) ) 
IF  (V)  120,110,110 

no  y  «  r  ♦  u 

120  h  •  y 

H  -  128  *  M  «•  16984 
IF  <M  -  Nl  140,140,130 
130  L  •  N  -  128*(N/128) 

1810  •  l  ♦  M 
CO  TO  150 
140  IB IC  -  N 
ISO  RETURN 
END 


FUNCTION  LLOWIKI 

LLOW  ■  K  -  12B*( K/1261 

RETURN 

END 


FUNCTION  TMOKIKi 
L  >  K  /  128  *  128 

TWQK  •  2-**L 
RETURN 
ENO 
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